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Abstract 

The Generalized Finite Element Method (GFEM) is a Partition of 
Unity Method (PUM), where the trial space of standard Finite Element 
Method (FEM) is augmented with non-polynomial shape functions with 
compact support. These shape functions, which are also known as the 
enrichments, mimic the local behavior of the unknown solution of the un- 
derlying variational problem. GFEM has been successfully used to solve a 
variety of problems with complicated features and microstructure. How- 
ever, the stiffness matrix of GFEM is badly conditioned (much worse com- 
pared to the standard FEM) and there could be a severe loss of accuracy 
in the computed solution of the associated linear system. In this paper, 
we address this issue and propose a modification of the GFEM, referred 
to as the Stable GFEM (SGFEM). We show that the conditioning of the 
stiffness matrix of SGFEM is not worse than that of the standard FEM. 
Moreover, SGFEM is very robust with respect to the parameters of the 
enrichments. We show these features of SGFEM on several examples. 

Keywords: Generalized finite element method (GFEM); partition of unity 
(PU); Extended Finite Element Method (XFEM); approximation; eondition 
number, loss of accuracy, linear system; Validation and Verification 

1 Introduction 

During the last decade, the Generalized Finite Element Method (GFEM) and 
the extended Finite Element Method (XFEM) - two approaches based on the 
Partition of Unity Method (PUM) - were developed independently and have 
been widely used to solve various types of problems. Only recently, it was 
clearly recognized that these two methods are same and were referred to as 
XFEM/GFEM ([13|). Hence we believe that it is interesting to briefly describe 
the early development of these methods. It was also recognized that, though 
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these methods have excellent convergence properties, the stiffness matrices asso- 
ciated with these methods could be ill-conditioned. In this paper, we especially 
address this issue and propose an easy modification, which we call the Stable 
Generalized Finite Element Method (SGFEM), that docs not have the above 
mentioned conditioning problem and is very robust. 

We start with a brief history of the early development of the methods based 
on PUM. Since this is history and not a survey, it is important to provide not 
only the publication date, but in addition, the submission dates for various 
papers, and pay careful attention to nomenclature for the various methods. 

Brief early history: The idea of adding non-polynomial basis functions into 
the trial space of the FEM started in 1970's ( [TOl [II [H] ) . However, these basis 
functions had global support and the associated stiffness and mass matrices lost 
their local structure. 

Three Special FEMs, which used non-polynomial shape functions, were pro- 
posed in [5] (1994, sub: Mar. 1992) to solve second order problems with rough 
coefficients. In particular, the shape functions used in the Special FEM #3 
have compact supports and are products of pieccwisc linear FE hat-functions 
and a non-polynomial function that mimic the special features of the unknown 
solution. This idea was further generalized with detailed mathematical theory 
and applications in the Ph.D. dissertation of J. M. Melenk 22](1995), where 
it was shown that the hat-functions could be replaced by any PU (with com- 
pact support). This method was referred to as PUM and PUFEM in f55]fl996. 
sub: Apr. 1996) and [n](1997, sub: Jul. 1995), respectively; these papers contain 
major results on the method and its application to the problems with highly 
oscillatory solutions, problems with solutions with boundary layers, differential 
equations with rough coefhcients, etc. 

The PUM was referred to as the GFEM in [47](2000, sub: Jul.1998), |48](2000, 
sub: Nov. 1998), [49] (2001, sub: Jul.2000), where the hat-functions were used 
as the PU (similar to the Special FEM #3 in 0). In these papers, GFEM is 
interpreted as an FEM augmented with non-polynomial shape functions with 
compact support, and it is shown that the use of only a few of these shape func- 
tions is enough to address the problems with singular solutions. Moreover, the 
idea of obtaining the non-polynomial shape functions by solving certain local 
problems is also introduced in these papers in the context of the analysis of a 
perforated plate. 

In a parallel development, but independent of [47], [48], [49], PUM with hat- 
functions serving as the PU was also investigated in [9](1999, sub: Jul.1998) 
and [35], (1999, sub: Feb. 1999) in the context of crack propagation problems. 
This method is similar to GFEM as it also uses the standard FE trial space aug- 
mented with non-polynomial shape functions. However, the major contributions 
of these papers is to show that the method does not need remeshing as the crack 
propagates. Also shape functions with jump discontinuities were used in these 
papers. The method was first referred to as the XFEM in the Ph.D. thesis of 
J. Dolbow [16](1999) and almost simuhaneously in [50](2000, sub: Scp.1999), 
[13] (2000, sub: Sep. 1999], and [15] (2000, sub: Sep. 1999). We mention that 
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XFEM was employed in [501 1131 [H] to address crack propagation problems in 
3-d, problems with branched cracks, and fracture in Reissner-Mindlin plates. 

Another idea similar to the PUM was used in the h-p Cloud method in 
[n](1996, sub: Jun.1995), [Tl](1996, sub: Apr.1996), where the shape functions 
were the products of a PU and polynomials. The goal of this method is to obtain 
h-p FEM like approximation without using a FE mesh, in the spirit of meshless 
methods. The use of the "customized function" (which mimicked the exact 
solution) for crack problems was also suggested in [35|(1997, sub: Dec. 1996), 
under this framework. Later, the hat-functions were also used as the PU in the 
h-p Cloud method in [3S](1998, sub: Dec. 1996). 

Lot of work has been done in the area of GFEM and XFEM since these early 
work, described above. We will comment on some of the recent developments 
near the end of this section. 

GFEM and the problem with conditioning: PUM is a flexible framework 
to design Galcrkin methods that accurately approximate solutions of variational 
problems. The framework involves (a) accurately approximating the solution, 
locally, using functions in a local approximation space, and (b) gluing the local 
approximations, using a PU, to construct a globally conforming approximate 
solution. The GFEM, which is a PUM with FE hat functions serving as the 
PU, retains the important flexibility of choosing the local approximation space. 
The efficiency of GFEM lies in the fact that is requires only modifying an 
existing FE code to incorporate special shape functions with compact support. 
The GFEM, with appropriate choice of special shape functions, leads to excellent 
convergence properties. However, the use of hat- functions as PU may result into 
almost linearly dependent shape functions in GFEM, and the stiffness matrix 
could be severely ill-conditioned; the ill-conditioning could be much worse than 
the conditioning of the stiffness matrix of the FEM. This results into the loss 
of accuracy in the solution of the linear system associated with the GFEM. In 
fact, the shape functions could be linearly dependent yielding a singular stiffness 
matrix. 

Various ad-hoc approaches have been developed in the literature to address 
this issue. For example, the stiffness matrix of GFEM was perturbed by an 
identity matrix of size e (small) in [47( I49| and an iterative method was used 
to solve the perturbed linear system. Preconditioning of the stiffness matrix, 
based on domain decomposition, have been recently suggested in |34| to address 
the conditioning problem. In |251 143] , a flat-hat PU (modified FE hat functions 
with flattened top) was used in the PUM instead of hat-functions. The use of 
flat-hat PU certainly avoids the problem of loss of accuracy in the linear system, 
but it requires developing a code from the scratch. 

Naturally, it is pertinent to ask if GFEM could be modified so that it retains 
the excellent convergence properties of the GFEM, and the loss of accuracy in 
the computed solution of the linear system of the modified GFEM is of the same 
order as that of the standard FEM. In this paper, we will show that the SGFEM 
has both of these features. We have chosen a 1-d problem to present the idea 
of the SGFEM primarily for the clarity of exposition and not to obscure the 
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analysis with details that are not directly related to the SGFEM. However, the 
ideas and the associated analysis (including the notational machinery) could 
be easily generalized to higher dimensions and will be reported in a future 
publication. 

Indicator of the loss of accuracy in computed solution of the linear 
system: Consider the linear system Ax = b, associated with FEM, GFEM, or 
SGFEM, where A is an n x n sparse symmetric positive definite matrix. Let 
X be the computed solution of the linear system, obtained from an elimination 
method encoded in a linear algebra package and the computations follow the 
IEEE standard for floating point arithmetic (with guard digits). Set ij :— \\x ~ 
x\\2/\\x\\2 - the relative error that measures the loss of accuracy in the computed 
solution. Tj depends on the round-off, but in general, it also depends on the 
elimination algorithm and its implementation in the package, the compiler, the 
processor, and the computing platform with single or multiple processors, rj is 
related to the relative error in the approximate solution due to round-off. 

We seek an indicator that reliably indicates the loss of accuracy in the com- 
puted solution, characterized by 77, and is practically independent of other fac- 
tors mentioned above. Let H = DAD, where Z? is a diagonal matrix with 
Da ~ . Define the scaled condition number A{A) of A by A{A) := K2{H), 

where K2{H) = Ibll-ff^^lU is the condition number of H based on the || • |j2 
vector norm. We hypothesize that A{A) is the indicator, which we formalize as 
follows: 

Hypothesis H: 

7? w Cn^A{A)e; /3 w 0, (LI) 
where e is the machine precision. 

We will elaborate on the precise meaning of the hypothesis and validate it 
in the Appendix, borrowing the ideas from the area of Validation and Ver- 
ification. The indicator A{A) will be used to compare various GFEMs with 
respect to the loss of accuracy in the computed solution, which will allow us 
to choose a preferable GFEM. In particular, we will show in this paper that 
^(^j^SGFEM-^ < j?(A'=^^^*0, where A^^'^^^*^ and A^''^^' are the stiffness matri- 
ces of SGFEM and GFEM, respectively, and therefore the SGFEM is preferable 
over the GFEM. 

Some current work in GFEM/XFEM: These methods have been used in a 
variety of applications. For example, XFEM has been used recently to address 
two-phase fluid flow problems (|19j). mechanical behavior of nano-structurcs 
(|20]). and heterogeneous material with random interfaces ([SS]); GFEM has 
been used to address heat transfer problems with sharp thermal gradient f|40]). 
grain boundary in polycrystals (53), and electromagnetic problems ([30]). Spe- 
cial shape functions for problems with locally periodic coefficients are con- 
structed in |31j that yield exponential order of convergence. Also local problems 
to compute the shape functions for problems with rough coefficients are con- 
structed in [3], and it has been proved that GFEM yields exponential order 
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of convergence. For an extensive collection of references in XFEM/GFEM, we 
refer to [23] . 

Organization of the paper: In Section [51 we give the model problem in 1-d. 
We intentionally chose the problem in 1-d so that we could communicate the 
main ideas of SGFEM, when applied to this problem, in a fairly general fashion, 
without the notational and other technical complexity associated with higher 
dimensions. We describe the PUM and GFEM, together with the approxima- 
tion results, in Section [3] and show the conditioning problem in GFEM on an 
example. In Section [U we first describe the SGFEM in a simpler setting, show 
that SGFEM retains the convergence properties of GFEM, and establish that 
the scaled condition numbers of the stiffness matrices of the SGFEM and FEM 
are of the same order. We chose the simpler setting primarily to communicate 
the main idea of the method and the associated analysis. We then describe the 
SGFEM and provide the analysis in full generality. We note that some of the 
ideas presented here could have been presented in a simpler fashion by using 1-d 
arguments. However we did not take this approach; the notations and frame- 
work of the analysis, developed in this section, could be easily generalized to 
higher dimensions. In Section [5J we applied SGFEM to three specific examples, 
namely, interface problems, problems with singular solutions, and problems with 
discontinuous solutions. In the Appendix, we discuss the validation of Hypoth- 
esis H and present many validation experiments. We note that the Appendix is 
a very important part of this paper 



2 Model problem 

Let il = (0, 1) and, for an integer k > 0, wc denote the standard Sobolev 
spaces by if'^(r2) with the norm || • ||//fc(n) and seminorm | • |i/<!(n); for fc = 0, 
7?°(r2) = L^(r2). We would also use the spaces H^{A), where A is a sub-domain 
of ri. Consider the variational problem 

ueH\n), B{u,v)=F{v), yveH\n), (2.1) 

where 

B{u,v):= / au'v'dx and F{v) -.^ / fvdx (2.2) 
Jn Jn 

such that F{1) = f dx = 0. We assume that the function a{x) is bounded, 
i.e., there are constants a, (3 such that 

< a < a{x) </3, V .t e O (2.3) 

We note that a{x) could be smooth, but it also could be rough. It is well known 
that the problem (|2.ip has a unique solution, up to an additive constant. 

We define the Energy norm, of u G H^{A), where A is a sub-domain 

of ri, by 

\\v\\'^i^^<^ -.^ Ba{v,v), where Ba{w, z) / aw'z'dx. 
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It is well known that the solution u of (|2.ip is also the solution of a boundary 
value problem (BVP), posed in the strong form as 

~[a{x)u']' f, au' (0) = au' (1) ^ (2.4) 

provided au' is differentiable. 

3 Generalized Finite Element Method (GFEM): 

Let iS be a finite dimensional subspaee of _ff^(f2). The Ritz-Galcrkin method to 
approximate the solution u of (|2.ip is given by 

Uh e S, B{uh, v) = F{v), VveS. (3.1) 

The solution Uh is unique up to an additive constant. We can obtain a unique 
solution by imposing a natural constraint on Uh, namely, Uh{0) — 0. 

A Partition of Unity method (PUM) is a Ritz-Galcrkin method, where iS 
is constructed employing a (a) Partition of Unity (PU) and (b) Local approxi- 
mating spaces. A Generalized Finite Element method (GFEM) is a PUM with 
special PU. We first briefly described the PUM. 

For a parameter h > 0, Let l'^ := {i € Z : < i < N}, where N = N{h) is 
an integer. For i e let := (a'\6'') C ft such that (i) Q = U^g/hW-', and 
(ii) any x G fl belongs to at most k of the open intervals lu^; k is independent of 
i, h. The open interval ut^ is called a patch. Subordinate to the cover {uj^}j^^jh, 
let {iVf lig^h be a C° PU satisfying 

J2 nH^) = 1, V X e r!, mth^^n) < C, diamK"}||(A^/')'|U^(n) < C, 

where C > is independent of i (for details, see [3 1331 HI)- 

On each patch , i G /'', we consider an (n^ + l)-dimensional space V/^ - 
the local approximating space, namely 

y> = span{^f' '}^4o, V^'" e H\u:.,) and ^^1''^ = 1, (3.2) 

where n^s are non-negative integers. The functions (p^J^'^, j > 0, arc carefully 
chosen such the functions in mimic the the exact solution u, locally in w^. 
We will further comment on this issue later. In the rest of the paper, we will 
write /, Wi, Ni, Vi, i/?^'' in place of I^, wf, N!^, V/^, f^j^''^, respectively, with an 
understanding that they depend on the parameter h. The PUM is precisely 
p.ip . with the finite dimensional space S is given by 

S^Yl = span{iV,; , < J < n„ z e /} := 5i + 52, (3.3) 

where 

rii 

= {c : c = E vo^^^m, 52 = {c : c - E E (3.4) 

iei iei j=i 
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and j/o', yj*' G M. The functions ip^j\ j > I, and the associated spaces Vi are 
sometimes referred to as enrichments and enrichment spaces respectively in the 
hteratm^e. We will refer to Si as the basic part of S; S2 will be referred to as 
the enrichment part of S. Moreover, we will refer to the Galerkin method with 
iS = iSi as the basic part of PUM. Thus every PUM has a basic part based only 
on the PU. 

We now present the main approximation result of PUM in the Energy norm 
(see [giMlli). 

Theorem 3.1 Suppose u S H^{il,). Suppose for i £ I, there exists € Vi and 
Ci > 0, independent of i, such that 

\\u - C1l2(„^) < Cidiam{uJi) \\u - CWs^lj,) and \\u - ChiiUi) < Ei- 

Then there exists w G 5 such that 

\\u~v\\,^n)<C[j:..^je^]'^\ (3.5) 

where the positive constant C depends on k, Ci,/3/a. 

It is immediate from Theorem 13.11 that the PUM solution Uh £ S — Si + S2 
of p.ip satisfies 

< mf <C[E.e/^f]'^': (3-6) 

where u is the solution of (|2.ip . It is clear from above that the global accuracy 
of the PUM solution Uh depends on how accurately the solution u of p.ip can 
be approximated by the functions in Vi, locally on the patches uji. 

We mention that in higher dimensions, the patches uji are subdomains, which 
can have quite general shape. Theorem 13.11 as presented above, is also true is 
higher dimensions. 

We now describe the GFEM. Recall that the choice of PU in PUM is arbi- 
trary. The GFEM is a PUM, where (a) the patches uJi are "EE stars" relative 
to a finite element (EE) triangulation of f2, and (b) the piecewise linear EE 
hat-functions Ni, associated with the vertices of EE triangulation, serve as the 
PU. 

Let N = 1/h and recalhng I = {i : < i < N}, let T := ^ ih : i e I. 
Let {Tk}k&i\{o} be the uniform mesh on £7, where := [xk-i,Xk] are the 

elements; Tk'.= {xk-i,Xk) is the interior of r^.. The points Xi are called the 
vertices of the mesh. The patches {wijig/ arc defined as uji :— (a;i_i, I'^+i), 
i = 1, 2, • • • , — 1; also luq := {xq, Xi) and ui^ := (xjv^i, 2:^). Eor i £ I, let Ni 
be the standard hat- functions associated with the vertex Xi ; the support of Ni 
is uJi. Note that luq = ti, lJ^ = and Ui ~ TiU Ti+i for i ~ 1, 2, . . . , — 1. 
ZUi is the EE star associated with the vertex Xi. Clearly, {Ni}i^j form a PU 
subordinate to the patches {uji}i^i. The associated GEEM is the Galerkin 
method (|3.ip with S = Si + S2 (see (|3.3[) ). Clearly iSi is the standard EE space 
of piecewise linear functions, and consequently, the basic part of the GFEM is 
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the standard finite element metliod (FEM) . Thus the trial space S of the GFEM 
is precisely the standard FE trial space, augmented with the space Thus 
GFEM could be implemented by incorporating enrichments into an existing FE 
code. The name GFEM was first used in [571 HH] to highlight exactly this point. 
The description of GFEM is exactly same in higher dimensions; it is based on 
the standard FE triangulation of fl. 

Remark 3.2 We note that we have considered a uniform mesh only for the 
simplicity of exposition; in fact, the ideas and theory in this paper could also 
be presented for locally quasi-uniform meshes, i.e., when Ci < |Ti;_|_i|/|rfc| < C2 

for /c = 1, • • • , — 1, with Ci, C2 > independent of k. ■ 

The accuracy of the GFEM (also PUM) solution depends on the choice of Vi, 

as mentioned before (see Theorcm l3.ip . The functions ip^^ G Vi (sec (|3.2p . (|3.3[) ) 
are carefully chosen based on the available information on the unknown solution 
u of (j2.ip to mimic the unknown solution locally in w^. Examples of Vi, suitable 
for specific applications are available in the literature (e.g., see [3]). We briefly 
mention some of the examples that we will consider in this paper: 

• If the unknown solution u is smooth in cj^, then the f^J^s are usually chosen to 
be polynomials in V^{uJi) and the associated spaces Vi are spaces of polynomials 
of degree n^. We note that rii could could be different for different i, based on 
the available information on u. 

• When a{x) is a piecewise smooth and discontinuous function (interface prob- 
lems), ip^^^s are chosen such that a[(^j'']' is smooth on cui. Clearly, (^'-'^ are 
continuous piecewise smooth functions with derivatives that are discontinuous 
at the discontinuities of a{x). 

• If the unknown solution u is singular, then ipj should be chosen as singular 
functions, mimicking the singularity of u. 

• If u is discontinuous at a; = c in the domain, then <Pj''s are chosen to be 
discontinuous functions on those w^s that contain x = c. We note however that 
problems with discontinuous solutions cannot be cast as (|2.ip : wc will address 
these problems in Section [5731 of this paper. 

Remark 3.3 GFEM provides a flexible framework to obtain various Galerkin 
methods. Many classical methods could be cast in this framework. For example, 
with rij = for z e / in GFEM (with S ^ Si) yields the classical FEM. 

Moreover, let s{x) be a function (could be singular) defined on fl. Consider 
Hi = 1 and ip^l^ = s{x)\^^ in the definition of S in ([53). Then GFEM, with 
Ui^ = 6 (a constant) for i G I, yields the classical "singular FEM" (see [JHl IHl 
[TTl [T^l Unmi] ) , where the standard finite element trial space is augmented by 
the global function s{x). Moreover, we note that rii in (|3.2p could be different 
for different values of i. In fact, one may use > only for a few patches coi, 
as needed for accuracy, based on the available information; for other patches. 
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rii = 0, i.e., Vi = span{l}. This idea was also discussed and implemented in the 
original GFEM papers [17 1 148 ] . 



3.1 Scaled condition number of the stiffness matrix of 
GFEM 

The stiffness matrix A of the GFEM is positive semi-definite. Even when the 
GFEM solution is naturally constrained with Uh{Q) = 0, i.e., when A is positive 
definite, the condition number K2(A) can be extremely large, specifically larger 
than the condition number of standard FE stiffness matrix, which is 0(/i~^) 
for second order problems. However, according to Hypothesis H, the scaled 
condition number .S(A) = K2{H) is a reliable indicator of the loss of accuracy 
in the computed solution of Ax = h. Recall H = DAD, where _D is a diagonal 
matrix with Da = A^- . We now present an example where .ft(A) is much 
larger than the scaled condition number of the standard FE stiffness matrix, 
which is again 0{h^^). 

Suppose a{x) = 1 in and let u = x" with 1/2 < a < 3/2, a ^ 1. Note 
that € H^{n) but a;" ^ H^{n). We consider a GFEM with n, = 1 and 
ifi^ :— i Cz I. From the definition of S, any u G 5 is of the form 

v{x)^ J2 a,N,{x)+"^hN,{x)x''; a„ 6, G M. (3.7) 

iG/\{0} 16/ 

We have set oq = to impose the constraint Uh{0) =0 on the GFEM solution 
Uh- It can be easily shown that ~ u, i.e. there is no approximation error. 

We let 77 := [ai, . . . , oat, &o, • ■ ■ e K^w+i, Then B{v,v) = At], 

where A is the (2A^+1) x (2A^+1) positive definite stiffness matrix of the GFEM. 
We note that Aii = |iVi|^i(„.) for 1 < i < iV and AAr+i+^.w+i+i = \N]X°'?hi(^.) 
for < j < N . Therefore by considering u G 5 of the form 

/ N N,(x) v^, NJx)x°' , ^ „x 

,67^0} fti \N^x"\HHn) 

it is easy to see that B{v,v) = rj^Hi], where H is as mentioned before. 

We consider a 17 G 5 of the form p.8|) with = for 1 < i < A^ — 1, on = 1, 
and 6i = for z G /. Then 

B{v,v) = J^v'^dx^l and ||7]||2 := |^ a^ + J2b^ = l. 
Therefore, 



ie/\{o} iei 

B{v,v) 



1 < Am, (3.9) 



where Am is the largest eigenvalue of H. 

Let g{x) G H^{Vl) be a non-decreasing function with g{x) = for < 
a; < 1/4 and < C < g{xi) < 1 for z > \N/2\. For h small enough, let 
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l/8<Xfe<l/4be the vertex closest to a; = 1/4. Clearly and gx" are in 
H'^{il), where := (1/8, 1). We now consider a w G 5 of the form (|3.8|) with 
tti = -g{xi)xf\Ni\Hi{Q) and bi = gixi)\NiX°'\Hi{n)- Then 

AT N 

v{x) = -^g(xOxfiV,(a;) +^5(a:.)iV,(xK. 

z— A; I— /c 

Thus i; = on [0, Xk]- Moreover, on r^, i > A: + 1, we have 

v\r^=-Iligx'^) + x'^Ilig), 

where is the linear interpolant of / on r^, interpolating at Xi^i and x^. 

Therefore, 

where we have used standard interpolation estimates. Thus recalling that w = 
on [0,Xfc], we have 

i=k+l 

l.2ri„„ct|2 I 11 Q!||2 ll„l|2 1 . /oi,2||| a\\\2 



Also, 



i=k i=\N/2'] 

where we have used that \N^x°'\jji^^^ > C/h for i > \N/2]. Thus using ([3TT0| . 
we have 

and hence, 

A™ <C/i4|||5a;"||p 
where Am is the smallest eigenvalue of H. Finally, from p.9p . we get 

Mly.^ Ml 

which is much bigger than the scaled condition number of the stiffness matrix of 
the standard FEM; we recall that the standard FEM is basic part of the GFEM. 
Thus from Hypothesis H, there will be severe loss of accuracy in the computed 
solution of Aa; = b. We will show this feature in the Appendix. 

It is interesting to note that using v £ S oi the form (|3.7p and following 
the same arguments as before, we can also show that the condition number 
K'2{A.) > C/i^^/|||gx"||p. We stated this property at the beginning of this 
subsection. 
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4 Stable Generalized Finite Element Method 
(SGFEM): 

A GFEM will be referred to as an SGFEM if the GFEM satisfies the following 
property: the scaled condition number ^(A) of the associated stiffness matrix A 
is of the same order with respect to h as of the stiffness matrix of the basic part 
of the GFEM. Since the basic part of any GFEM is the standard FEM, therefore 
a GFEM is an SGFEM provided J^(A) ~ 0(h~'^) for second order problems. 
As mentioned before, we will present the analysis for uniform meshes. However, 
the analysis is valid for locally quasi-uniform meshes. 

We first present a particular example highlighting the ideas and results re- 
lated to SGFEM in a simpler setting. 

4.1 An example of the SGFEM: 

Let a{x) = 1 in (|2.ip and suppose the solution of (j2.ip is smooth, in particular 
let u € H^{Cl). Since the solution is unique up to an additive constant, we seek 
u with m(0) = 0. It is well known that a function in H'^{il) could be accurately 
approximated, locally in uji, by polynomials of degree 2; recall that the patches 
Wj have been defined in Section [31 Based on this information, we consider 
Vi = span{(^'*'}|^Q (i.e., rii = 2), where tp^^'^ = {x — Xi) and ip2^ = {x — x^)^, for 

0<i<N. RccaU that (/^[j'' 1. Thus = V'^iuJi). 
We let 

T^:=vf-I..{^), where2;..(^W):= ^ ^{x,)N,l^; 

l-l<k<i+l 

{'pf' ) is the piecewise linear interpolant of (/s'*' on the patch uji. We adjust the 
operators X^^g andX^^^ ; they interpolate at {xq, xi \ and {xn-i,xn} respectively. 
We define a modified local approximation space Vi = span{^'*'}|^Q, associated 

with Vi. Clearly, = for j = 0, 1 and thus Vi = span{^2'}- 

It is well known (see [33l|49]) that the scaled condition number of the stiff- 
ness matrix of the GFEM, with Vi as the local approximation spaces, could be 
extremely large or even unbounded. We will use the GFEM with Vi precisely to 
address this issue, and show that the GFEM based on the approximation space 

S = S1+S2, with 5i = ^ OiN, and ^2 ^^NiVi 

tei\{Q} iei 

is an SGFEM. Note that w(0) ~ for all v <E S. We have chosen ap = in the 
definition of Si to impose the constraint w/i(0) = to obtain a unique GFEM 
solution Uh e S. 

It is easy to check that the assumptions in Theorem 13.11 hold: in fact, there 
exists e Vi such that ||u — ^ Ch'^\u\}j3i^^.y Therefore it is clear from 

(|3.6p that ||u — Uh\\£{n) = C>{h'^), where Uh is the GFEM solution, based on 
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S = Si+S2 (recall ^2 = Y.tei ^^^O- We first show that the GFEM based on 
S ^ Si + S2 also yields the same optimal order of convergence. 

Proposition 4.1 There exists a v E S = Si + S2, such that 

\\U - < C/l^|w|i/3(n), 

where the positive constant C independent of h. 

Proof: Since u E H^{u!i) ioi < i < N, it is well known that there exists 
f eV,^ V^iiOi) such that 

\\u-e\\£(u.)<ChMHHu.,)- (4.1) 

Let Z/iU = X^ie/ u{xi)Ni. It is clear that IhU = 2i^.u on uJi, Therefore using 
standard interpolation results, we have 

ii(u-x„^.)-(r-i.,r)iiL2K) = \\{u-e)-:c.Au-e)\\LH.,) 

< CdmiR{uji)\\{u- C)\\£{u,), 

and similarly, 

\\{u-I,,u)-{C-I^A')\\£iuO < C\\u-C\\£iu.^)<Ch^\u\H^u.^y 

Let w := u—Ihu; clearly w G H^{rt). From above, X^^^^* G Vi approximates 
w locally in w^. Therefore, from the Theorem 13.11 there isvGS2 such that 

Let V = I/jU — u{xq) + V. Since {Ni}i^i is a PU, we have I^u — u{xo) = 
Z]i6/\{o}["(^«) ~ 'U'{xo)]Ni <E Si. Thus V e S and using (|4!2)) . we get 

ll"- w||£(f2) = llw' - t'lUcn) < C'/iVlff3(a), 

which is the desired result. □ 

Using Proposition l4.1l we immediately get that — u/iH^jo) = 0{h'^), where 
Uh is the GFEM solution based on iS = iSi + 1S2. We also note that we approx- 
imated u —IhU by the functions in Vi in the proof of Proposition l4.1| - this is 
the main idea of SGFEM. Later, we will further comment on this issue. 

We now address the scaled condition number of the stiffness matrix A asso- 
ciated with the GFEM based on S = Si + S2- With a suitable ordering of the 
shape function of S, the matrix A is of the form 



All A12 
A21 A22 



(4.3) 



where A^j are block matrices. The matrix An = {^(-^i, -^j)}ijg7\{o}: which 
is the stiffness matrix of the basic part of GFEM, is the standard N x N EE 
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stiffness matrix. The {N + 1) x {N + 1) matrix A22 is of the form A22 = 
{B{Ni(p2\ Nj'ip2^)}i,j(zj. Also A21 = A|2- For the clarity of notation, we 
win write A22 = {(^22)^ }.fj=iJ where M = N + 1. Note that (A22)jj are 
associated with the vertices Xj-i, respectively, and the GFEM solution Uh is 
computed by postprocessing. We remark that, in general, M will vary based on 
the application and {A22) jj will be associated with some vertex x^q). 

We first note that Tp^^Xj) = Oforj = i — + Therefore it is easy to 
show that Si and 52 arc orthogonal in the inner product B{-, •), i.e., 

B{vi,V2) = 0, y VI eSi, V2 eS2. (4.4) 

Thus it is immediate that A12 and A21 in (j4.3|) are "zero-matrices" . 

The matrix An is tridiagonal and is constructed by the assembly process 
from the element stiffness matrices A^^', for the element — [x^^i, Xk], k = 



1,2, 



A(k) 

, N. the matrices Al^' are given by 



A 



(fc) 
11 



^11 , ^11 



h 



1 -1 
-1 1 



, 2 < /c < TV, and A'{^ := [I]. (4.5) 



Similarly, the matrix A22; which is also tridiagonal, is constructed by the as- 
sembly process from the element matrices 



A 



(fe) _ 



BrANk-i^[' '\Nk^i^[^ '1) Br,{NwP,Nk-i^^^ '1) 



BrdNk-l^^t'\Nkip'2^') 



(4.6) 



for the element r^, k = 1,2,- •• ,N, and Br^{w,v) := J^^au'v' dx. A direct 
computation yields 



^22 ~ ^22 , ^22 



2 


1 


15 


30 


1 


2 


30 


15 



It is easy to check that the matrix A'"^ is positive definite (the eigenvalues are 
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and i) and thus 



Il2/||^ Vy = (yi,y2)e 



10""" - " " - 6 
We now consider the diagonal matrix D = rfia(7(Di, D2) with 



Di = dmg(di), di = m~^^^ {2-^^^ , ■ ■ 
D2 = diag{d2), da - m~^/^(l, 2-^/2, 

where mi = I//1 and m2 — 2h^/15. 
Wc next define 



, 2^1/2, 1)^GM^, 
•• ,2-1/2,1)^ eM"+i, 



DAD 



DiAuDi 
D2A22D2 



An 






A22 



(4.7) 



13 



where An = DiAnDi and A22 = D2A22D2. Clearly An and A22 are N x N 
and (iV+ 1) X (7V+ 1) tri-diagonal matrices, respectively. The diagonal elements 
of All and A22 are equal to 1. Consequently, diagonal elements of A are equal 
to 1 and the scaled condition number ^{A) of A is K2{A). 

Proposition 4.2 Suppose .S(A) be the scaled condition number of A and let 
Amin(Aii), Xmax{Aii) be the smallest and largest eigenvalue of An, respec- 
tively. Then 

< m) < ^(Aii) ".^f;'^f^^"'-f"| , 
min{l,6i/Am„(Aii)} 

where Ci =3/4 and C2 =5/4. 

Proof: Let z = (z^za)^ G M^w+i^ where Zi G and Z2 € M^+i. Then 

Z^AZ = (DiZif Aii(DiZi) + (D2Z2f A22(D2Z2) 

= zf AiiZi + Z^A22Z2. (4.8) 

Let Z2 = (2/1,2/2, ■ ■ ■ ,2/JV+i)^, then D2Z2 = ^ (2/1, 2"^2/2, ■ ■ ■ , 
2"2yjv,2/Af+i)^, where m2 = 2/1^/15. We define Z2,fe := {yk,yk+iV, and 

Z2,i := 7712^^^^(2/1, 2^^/22/2)'^, Z2.W := 77i2"^^^(2"^/^2/w, 2/Af+i)'^, 
Z2,fe := (27772)~'/'(2/fe, 2/fc+i)^, for fc = 2, • • • , - L 

Recalling that A22 could be obtained from the element matrices through 
the assembly process, using (|4.7p we get, 

AT 

Z^A22Z2 = (D2Z2)^A22(D2Z2) = ^Z2,fc^A^2^Z2,fc 

fc=l 

u3 ,3 N+1 

< yEil^-^ll^ = yE-.-V = ^l|z2f. 

k=l 1=1 

Similarly from (j4.7p . wc also get 

3 

-IIZ2IP < Z^A22Z2, 

and therefore from (|4.8p . 

zf Aiizi + C1IIZ2IP < z^Az < zf Auzi +C2i|z2|p. (4.9) 

It is clear from above that 

z^Az > zf Aiizi +Ci|lz2f 

(Aii)i|zi||2 + Ci||z2jp 
> min{Ci, A™„(Aii)}||z|p, 
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where Ci := |. Therefore, 

)} = Ami„(Aii)min{l,Ci/Amm(Aii)}. 
Similarly from the upper bound of (|4.9p . we can show that 

AmQx(A) < Xmax{Ml)ma.x{l,C2/Xmax{Ml)}, 

where C2 = |. Thus 

A{A) - K2(A) - "^'""^(^^ < ^ma:z(Aii) maxjl, C2 / Xmax{Ml)} 

A„„i(A) Amm(Aii) min{l, Ci/Ami„(Aii)} 

max{l,C2/Amax(Aii)} 



i^(An)- 



i{l,Ci/A„,,„(Aii)} 



where ^(An) ~ K2(Aii). Thus we have the required upper bound of ^(A). 

Now let Zi be an eigenvector of An associated with A„ia2-(Aii). Also let 
Z2 = 0. Then from (|4.8[) . we have 



z^Az = A„ax(Aii) llzif = A„a^(An) \\zf, 

and therefore, 

Amaa:;(A) ^ Aj7^aa;(Aii). 

Similarly, considering zi to be an eigenvector of An associated with Amm(Aii) 
and Z2 = 0, we have 

z^Az = A„i„(An) !|zi|p = A„„j(Aii) \\z\\^, 

and therefore, 

AmMi(A) < A,„i„(Aii). 

Now, 

ATrjin(A) A^i2^i(Aii) 

which is the required lower bound of ^(A). □ 

The Proposition 14.21 establishes that M.{A) sa R{Aii), i.e., the scaled con- 
dition numbers of the stiffness matrices for the GFEM and the basic part of 
the GFEM are of same order. Thus the GFEM with S = Si + S2 is indeed an 
SGFEM. 

Remark 4.3 We note that the orthogonality of the spaces Si and ^2 was essen- 
tial in proving Proposition l4.2l This property does not hold in general. Later we 
will define a notion of "almost orthogonality" of Si and 52 , which will address 
this issue. ■ 
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Remark 4.4 The inequahty (|4.7p played an important role in obtaining Propo- 
sition |4]2l This property depends on the functions in Vi. For general approxi- 
mation spaces Vi, we need an assumption that will be presented later. ■ 

Remark 4.5 SGFEM uses Vi as the enrichment space, which is a modification 
of Vi. Other modifications of V have been reported in different contexts. For 
example, the shifting modification, namely, Tp^^\x) = (p^j\x) — ip^j\xi), j > 0, is 
used in XFEM in the context of approximation error as well as enforcement the 
Kronecker delta property (see [23j).B 

4.2 SGFEM and its analysis: 

We now present the SGFEM for (gH]), with a e L°°{n) and < a < a{x) < p. 
Moreover, suppose it is a priori known that u(0) = (we will further comment on 
a priori information later). We consider the uniform mesh {Tk\k£i\{ii} with the 
set of vertices T, as described in Section [31 Recall that the hat function Ni and 
the patch uji is associated with each Xi S T. We will refer to {x^-i, x^, ajj+i} as 
the vertices of oj^; the vertices of cjq, are {a;o,cci}, {xn~i,xn}, respectively. 
Let 

Ti, T2 C T; Ci card(ri), C2 caix^Ta); Ci, C2 < + 1- 

We define iSi = X^a; sTi '^i-^i^ e M; 7i will be referred to as Si-relevant set of 
vertices. We consider 7i = {xi e T : 1 < i < A'^} as in the example in Section 
14.11 For other choices of 7i , we refer to Remark 14.61 

For Xi g T, let Vi = span{ }"1q C H^{ui) such that there exists e Vi 
satisfying ||u - C'll^j"^) < for all C Ui. Clearly, \\u - C\\£(LOi) < 2ei. We 
consider the modified space Vi ~ span{^^''}"4i7 where 

I<^.u is the piecewise linear interpolant of v £ H^iuji) on the patch uji based on 
the vertices of w^; we adjust X;^^ andX^^^ accordingly as before. It is important to 
note that if for some Xi G T, = G H^{uji) : ^\r^ G V^irk) for aU C uJi}, 
then Vt = {0}. Also Ti^k) = with fc = « - 1, i, i + 1 for aU J' G Vi. We 
refer to a patch uJi as enriched if Vi 7^ {0}. Let T2 := {xi G T : is enriched} 
and define ^2 — eT2 ^^^^ referred to as the S2-relevant set 

of vertices. In Section 14.11 we chose T2 ~ T. We will present examples with 
(2 << N + 1 (i.e., only few patches enriched) later in the paper. 

Remark 4.6 The sets 71, 72 C T provide a framework to address numerical 
treatment of many applications. Selection of both sets depends on a priori 
information on the problem and its solution. Selection of 72 will be apparent 
from the examples in Section[5l Suppose To d T contains all the vertices Xj G T, 
where it is known a priori that u{xj) — 0. We choose 7i = T\To. Typically, 
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Ti will not contain any boundary vertex with homogeneous Dirichlet condition. 
However 7i may exclude other vertices in T based on a priori information. For 
example, let f{x) = X^felo "-fc cos[27r(2fc + l)a;], a{x) = 1, and suppose it it known 
that u(0) = 0. Then m(1/4) = m(3/4) = 0, and the vertices Xj ^ 71 if Xj = 1/4 
or Xj = 3/4. Thus we can accommodate many a priori information in this 

framework. Only for simplicity, we have considered To = {a;o} in this section. ■ 

Wc now consider a GFEM with 

S^Si+S2= a,N, + N,V^ . (4.10) 

Note that v{Q) = for aU v e S. We will show that this GFEM is an SGFEM, 
under certain assumptions on the space 52, which we will present later. We 
mention that Ti and T2 are called Si and ^2 relevant vertices, respectively, 
since the degrees of freedom associated only with these vertices appear in the 
GFEM. _ 
Wc first present an approximation result for the GFEM with S = Si + S2- 

Theorem 4.7 Let u G H^{n) be the solution of (|2.ip . Suppose for each Xi € 
72, there exists ^* e Vi and Ci > 0, independent of i, such that 

\\u - I^^u - rilL2(^^) < Cidiam{LOi)\\u - I^.u - ril£(t^i): 
and \\u — Icji" — — ^i- Then there exists v G S ~ Si + S2 such that 

\\n~vhin)<C{ J2 Y ''f- (4.11) 

XieT\T2 X,£T-2 

Proof: Let XhU = '^j..^q-u{xi)Ni be the piecewise linear interpolant of u. 
We note that ThU = Ii^.u on cji. Define w := u—ThU and let v := 6T2 -^j? ^ 
^2. Then recalling that {Ni}xi£T is a PU, we have 

Therefore 

\\y^-nlia)<C[\\ J2 + E ^'(--^"^)C(a)- (4-12) 

X,£T\T2 x,eT2 

We first address the last term of (|4.12p . Using the fact that x G is in at 
most two patches uji, Wi+i, we see that the sum '^x-eT^^^i^ ~ ^')]' 
most two terms for any x £ fl. Using this observation, the assumption that 
ll^illL'°(f2) < C[diam{wj}]~^, and the hypothesis of the Theorem, we can show 
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that 



EAT/ pi\ ^ I'™ ^'Ili^ctJi) 



< E ii«^-eii^(..)< E (4-13) 

(We refer to the proof of Theorem 3.2 in [4] for details of the argument lead- 
ing to ()4.13p ). Using exactly same argument and the interpolation estimate 

lkllL2(a;,) = \\u-TuiU\\L^(uj,) < Ch\\u-I^^u\\£(^^.-), we get 

II E E !l«-^-."ll£(c..)- 

xieT\T2 xieT\T2 

Therefore, from (|4.12p and (|4.13p . we have 

\\^-niin)<c[ E E ^?]- 

Xi€T\T2 Xi^T2 

Finally, writing w = u —T^u and setting u = T;,m + U e 5i + we get the 
desired result. □ 

We mention that unlike in Theorem 14.11 we did not assume 72 = T in The- 
orem 14.71 We further note that XhU for u E {Vl) is not defined in higher 
dimensions, since the point values of u, in general, do not exist in higher dimen- 
sions (in contrast to 1-d). However, using a generalized interpolant based on 
the average of u in a ball around the vertices Xi , the proof of the above result 
can be easily generalized to higher dimensions. 



Remark 4.8 From the proof of Proposition 14. 1[ it is clear that accurate local 
approximation of u — ThU by functions in Vi is crucial to obtain the desired 
result. This is the main idea of SGFEM - the spaces Vi are constructed such 
that the functions in Vi accurately approximate u—ThU in uji. This is in contrast 
to the standard GFEM, where the functions in local approximating spaces Vi 
accurately approximate it in w^. ■ 

Remark 4.9 We note that V^ = {0} for Xi € T\T2- If u G H^{n) is locally 
smooth, namely, u G H'^{u!i) for Xi G T\72, then ||m — < Ch\u\H-^(^^.-) 

for Xi E T\72 and (|4.1ip could be written as 

\\u-v\\,^n)<C{h^ E \^\W.}+Y.''V^'- (4.14) 

XiGT\T2 Xi£T2 

By incorporating the available information on the solution u in V^, for Xi e 72, 
we can have = 0{h), and consequently, \\u — u||£(r2) — 0(h). The set 72 can 
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be chosen adaptively with respect to a prescribed tolerance, which we do not 
elaborate in this paper. ■ 



Remark 4.10 A rate of convergence of 0{h) for various problems have been 
reported for the Corrected XFEM (which is also a GFEM); see e.g., [22]. How- 
ever, for the crack propagation problems, the enrichment spaces Vi in XFEM 
requires the use of a ramp-function to obtain the 0(h) rate of convergence. 
In contrast, the GFEM based on 5 — 61+82 docs not require the use of a 
ramp-function to obtain the rate of convergence of 0(h). 

We now address the scaled condition number of the stiffness matrix of the 
GFEM with S = Si + S2. For clarity of the exposition, we will present the 
analysis for the case when rii = 1 i.e., Vi = span{^^*'} . The analysis for 
general Ui is similar. 

As in the example presented in Section 14. 1[ the stiffness matrix A is of 
A-ii A.\2 

the form A = . . , where An = {B{Nt, Nj)}x^,x^eTi is the (1 x Ci 
A21 A22 J 

stiffness matrix of the basic part of the GFEM. Let Di be a diagonal matrix 

— 1/2 

with (Di)ii = (Aii),jj . Glcarly, the diagonal elements of 



An :=DiAnDi 



(4.15) 



are equal to 1. 

The matrix A22 plays a central role in our analysis and depends on elements 
that have been enriched. We will refer to an element Tk = [xk-i, Xk] as enriched 



if (a) Xk-i e T2 and ip^^ ^^\r^ 



^ 0, or (b) Xk &T2 and ip\ 
:= {tu : Tfc is enriched} 



^ 0. Let 



The matrix A22 is constructed by the assembly process using the element stiff- 
ness matrices defined only on g JCem - 

We now address the structure of the element matrices A^22 i'^ detail and set 



up some notions and notations that will be used in the analysis. We denote 



the vertices of the element as x^ 



(fe) 



Tk S ICenr- The element stiffness matrix A22 is of the form 



Xk-l 

(k) 



and 



Xk] we consider only 



A. 



(k) _ 
22 ~ 



"11 "12 



-"12 



-"22 



(4.16) 



where = B., (A^fc-2 



.lk-2+i 



,Nk-2+](P 



1 < i,j < 2. 



A^2 is associated with the vertices x\^'' — x^-i and x^ 



lib\1',bi2^ > 0, then ^2' 

is 2 X 2 and we say that the local stiffness matrix 



(fc) 



a diagonal matrix D^'^' = diag{(5j , } with \5. 
diagonal elements of 



2 
(fe) 
2 



Xk- We define 
> 0, such that the 



^22 
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are of equal to 1 or 0(1), independent of h. Clearly, (5^'"\ (Jj'^'' are associated 
with vertices Xk-i,Xk respectively. 

On the other hand, if 622'' = (i.e., Ipf^ \r^=Q and consequently 6^2'' ~ ^21 — 
0) in (|4.16p . then the local stiffness matrix ^22'' = I^ii''] of size 1x1 and is 
associated only with the vertex x^*^^ = Xk-i- We define D^'''> = [s[''^], where 
jj'^-' = {b^^^}~^^^; jj'^'' is associated with the vertex x^'^-' = Xk-i- Similarly, if 
^ri = in (|4.16p . then the local stiffness matrix ^22'' = [^22''] associated only 
with the vertex x^''^ = Xk- Also [(jf''] with jf'' {bil^}~^^^ associated 

with the vertex x^^ ~ Xk- Let ^'^'^^ be the number of vertices associated with 

the local stiffness matrix ^22'' • Thus the size of ^22'' ^ note that q^^^ 

is either 1 or 2 with our assumption = 1. 

Recall that A22 is obtained by the assembly process using the element stiff- 
ness matrices ^22^ the size of A22 is C2 x C2- Let c = (ci,C2,-- - ,2^2)5 then 

c^A22C= M'^]''42^^^'^ (4-17) 

where c*^*^-' G M'"''"' . Moreover, the components of c^'""^ are also the components 
of c that correspond to those vertices of tu that are associated with ^22'' ■ For 
example, if ^n', 622'' > in ^22'! then as mentioned before, the vertices x^^'^ = 
Xk-i, x^2^ = Xk are associated with ^22"*. Suppose the components Cj(fc)_i, Cj{^k) 
of c are associated with the vertices Xk-i^xu, respectively, of r^. Then c*-*^-* = 
[cj(fe)_i, Cj(fc)]^. Similarly, if ^422'' = [^n''], then ^22'' associated with x'^i^ and 
(jC^) = [cj(i.)_i] - a vector with one component. Later in our analysis, we will 
use (|4.17p with a particular vector c and S-'''' as defined above. 

Next we note that each vertex Xi of the FE mesh is associated with a FE 
star - union of all elements Tk C w,; (equivalently, union of all elements Tk with 
common vertex Xi). For G 72, we define Ki := {rk G /Cenr : C Wi}. For 
e 72 and Tk € /Ci, we set the index 1 < fc) < 2 as follows. We first note 
that fee + 1}. For k ^ i, we set l{i, k) = /(i, i) = 2 and for fc = i + 1, we 
set /(i, fc) = i + 1) = 1. Thus fc) is the index such that x^^^j^-^ = Xi] note 
x'l^^ i^s^ may not be associated with ^22''. We define 

JC* := {Tk € JCi : a;||^^j,j is associated with A22'}. 

Thus K.* is the set of Tk E K-i such that f^i'lrt ^ 0. For G 72, we define 

A. := E ['5^S.)r^ (4.18) 

which will be used later in our analysis. 

Each diagonal element of A22 is associated with a vertex in 72. Let {A22) j^j^ 
be associated with G 72. Moreover, we note that {A22) jij^ — J^TkeK* ^1*^^ 



l(i,k),l(i,k)^ 



20 



where 6pq^ was defined in (|4.16|) . Tlius (A22)jiji > for all Xi € T2 (i.e., all 

the diagonal elements of A22 are positive). We now define the diagonal ma- 

2/2 

trix D2 = diag{di, (^2, • • • ,^^3} with dj = {A22)jj , 1 < j < C2- Note that 
dj- — {A22)j.j. is associated with a;^ G 72- Clearly, the diagonal elements of 

A22 := D2A22D2 (4.19) 

are equal to 1. Define the diagonal matrix D :— diag{Di,D2}. Since the 
diagonal elements of An, A22 (see (|4.15p . ()4.19|) ) are equal to 1, the diagonal 
elements of ^ ^ ^ _ 

" All A12 



A := DAD = 



A21 A 



22 



(4.20) 



are also equal to 1. Also A12 = D1A12D2 and A21 = Ai2- 

We will show that the GFEM with S ^ Si +S2 is an SGFEM, under the 
following assumptions on the local approximation spaces Vi and the enrichment 
part of S, namely iS'2. 

Assumption 1 The spaces Si and S2 are almost orthogonal with respect to the 
inner product B{-.-), i.e., there exist constants < Li,C/i < 00, independent of 
h, such that 

Li{\\vi\\j(^n) + \\v2\\l^n)} < l^ivi +V2,vi + V2)\ < + \\v2\\j(^n)}, 

for all vi G Si and V2 G ^2. 

Assumption 2 For Tk G ICenr, there exist constants < L2, U2 < 00, indepen- 
dent of k and h such that 

L2\\[D^'^]''^f < x^4'2^x < U2\\[D^'Y'^f, V X G R^''\ 
where the diagonal matrices Z?'*^' have been defined before. 

Assumption 3 For xi G T2, there exist constants^ < L^^Uy, < 00, independent 
of i and h such that 

L3 < (A22)-J,A, < [/3, 

where {A22) jij^ is the diagonal element 0/ A22 associated with Xi, and A; is as 
defined in (|4.18p . 

The following result is an easy consequence of Assumption [TJ 

Lemma 4.11 Let x = {(^ ^rf Y G where i G K^^ and 77 G M'^^. Then 

there exist positive constants Li and Ui, independent of h, such that 

Li [(^Aiii + ?7^A22?7] < a^^Ax < Ui [C^Anf + ?f A22ry] , 

where A, An and A22 are matrices defined before. 
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Proof: Let ^ = {$,i)x,eTi and 77 = {r]i)x^fzr2- Consider vi = Y^x^en^i^i ^ '^1 
and V2 = Z]a;ier2 '7»^#i' "= '^z- Then + V2,vi + V2) = x^Ax, B{vi,vi) = 
^"^Aii^, and B{v2,V2) = ri^A22'rj. The desired result is now immediate from 
Assumption [TJ □ 

Theorem 4.12 Suppose the Assumptions^^ [H and\^ are satisfied. Let A he 
the stiffness matrix of the GFEM with S = Si + 82- Then 

^ KiMi) < ^(A) < KiMi) El ^^ax {l,^.C.3/A„. (An)} 

Ul Li mm|l,L2i3/Am»n(Aii)| 

where Amm(Aii), Amaz:(A) are the smallest and largest eigenvalues, respectively, 
of the matrix An defined before. 

Remark 4.13 This result shows that under the Assumptions [U [21 and[3l the 
scaled condition numbers of the stiffness matrices of the GFEM with S = iSi+iS2 
and the basic part of the GFEM arc of the same order. Thus the GFEM with 
5 = 5i + 52 is indeed an SGFEM. 

Proof: Let z = (zi,Z2)^ G M^l+'^^ where Zi e R'^^ and za G R'^". Then from 
the definition of A (see ^^), we have z^Az = z^DADz = (Dz)^A(Dz), 
and since Dz = [(DiZi)-^, (DaZa)"^]"^, from Lemma [4. Ill we get 

Li [(DiZi)^An(DiZi) + (D2Z2)^A22(D2Z2)] < z^Az 

< Ui [(DiZi)^An(DiZi) + (D2Z2)^A22(D2Z2)] . (4.21) 

Let Z2 = (/i, /2, ■ ■ ■ , /c2)"^ and consider D2 = diag((ii, ^2, ■ • • , dc-i) with di = 
{^i-iUi^"^ as defined before. Then D2Z2 = (c?i/i, ^2/2, • • • id^^f^^)'-'^. Recah 
that dj. is associated with G 72- Consequently, dj-fj^ is associated with 
Xi G T2. 

Consider an element G /Cenr ■ Following the notation given after (|4.17p , let 
^2^^ ■= (D2Z2)*-'^'' G R""' ' such that the components of Zj'^'' are the components 
of D2Z2 corresponding to the vertices of Tfc associated with ^22^ Now from 
(|4.17p and using Assumption [2J we have 

(D2Z2)^A22(D2Z2) = ^ zf A^'^^z^ > L2 J] \\[D^''^]-'z!^r . (4.22) 
We note that if D^*^) = diag{(5f ^''^^l. then 

where Z2 — [djfk)-i^ fj(k)-i]'^ following the notation given after (|4.17p . Simi- 



-''j(^k) — l 1 J j(k)- 

larly, if i^C^) = [s[% then || [i^('=)]-iz^||2 = [S[''Y^[d,^k)-i?[.fm-i? ^ ^nd if 
DC^) = [Si% then i|[i?('=)]-iz^-f = [Si'YHdm^ifm]' 
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Now, it is important to note that if Ji := {dj(^k)-i.fj{k)-i,dji^k)fj(^k)}rkeK,„r 
(where the repeated elements appear only once) and '■= {djifji}xieT2J then 
Jt^ = J2. Thus from wc have 

(D2Z2)^A22(D2Z2) > ^2 ^»4/l ^ ^2^3112211', (4.23) 

where we used Assumption [3] to get the last inequality. Similarly, we can show 
that 

(D2Z2)^A22(D2Z2) < C/2 C/3 1 1 Z2 1 1 ^ 

Therefore from (|4.2ip and using the definition of An, we get 
Li zf Aiizi +L2L3IIZ2IP < z^Az<C/i zfAiiZi + L/2?73||z2||^ . (4.24) 
Now from the lower bound of z'^Az in (|4.24p . we have 
z^Az > Li A™„(Aii)||zi|p + L2i3||z2i 
and therefore, 

A™„(A) > ii A™„(Aii) min {l,i2i3/Amm(Aii)}. (4.25) 
Similarly, using the upper bound of z^Az in (|4.24p . we can show 

Amax(A) < UiXmaxiMi) max { 1 , [/j [/a/ A,„a2; ( An ) } . (4.26) 

Thus from (|4?25l) and (|4?26)) . we have 

^(A) ^ WA) ^ ^^^^^^ ^ max{l,t/2t/3/A„..(An)} ^^^^^^ 
A„„„(A) Li mm{l,L2i3/Amm(An)} 

which the required upper bound. The required lower bound could be obtained 
by following the exact arguments in Proposition 14.21 and (j4.24p . Thus we get 
the desired result. □ 

We mention that the notions and notations developed leading to Theorem 
I4.12l can also be extended to higher dimensions. An element will have He vertices, 
e.g., Ue could be 3 or 4 in 2-d. And the element stiffness matrices could be 
at most rig X rig. The assembly argument (|4.17p could be easily generalized to 
higher dimensions. For a given vertex Xi and an enriched element Tk in the FE 
star associated with Xi, the index l{i,k) will again represent the local index of the 
vertex a;||^^''j,^ of that coincides with Xi, i.e., x^i^^ = Xi. The expressions for 
Ai, (A22)m and the Assumpsions [TJ [21 [31 are exactly same in higher dimensions. 
Using these notions, the proof of Theorem 14 . 1 2 1 can be easily extended to higher 
dimensions. The approach presented here can also be extended for elasticity 
equations etc. We note however, the notations become a little more involved if 

Tl i > 1. 
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Remark 4.14 We now make comments on the assumptions. The Assumption 
[T]is always satisfied in 1-d. Let Bq{u,v) f^^u' v' dx. Since Tp^^^Xk) = for 
k = i — + 1, it can be easily shown that i?o('*^ii ^^2) = for all vi G Si and 
V2 £ S2- Therefore, 

B{vi + V2,vi + V2) > aBo{vi + V2,vi + V2) 

= a[Boivi,vi) + Boiv2,V2)] > + ^2||^(o)]- 

Similarly, we can show that 

B{vi +V2,vi +V2) < ^iWviWl^n} + h^Wl^n)]^ 

and thus Assumption [T] is satisfied with Li = ^ and L2 = ^- In higher dimen- 
sions, this assumption has to be checked. 

Assumption [5] is equivalent to i2||y|p < y'^A'^^y < U2\\y\\^ for all y G M^""' . 
Thus A22 is uniformly positive definite in k and its eigenvalues are uniformly 
bounded. 

It is always possible to choose the diagonal matrix Z?''^^ such that Assump- 
tion [3] is satisfied. For example, it is easy to check that Assumption |3] is satisfied 
with L3 = [/g = 1 by choosing = diag{5^ S^''^} with ^ = )"^/^ The 
Assumption [3] is trivially satisfied with L3 = U3 = 1 when D^'^^ is a 1 x 1 matrix. 



Remark 4.15 As shown in the Appendix, the implementation of the SGFEM 
does not require scaling the stiffness matrix, i.e., the linear system involving 
the stiffness matrix A, and not scaled version A, is solved. The scaling was 
used only to define J^(A) and to study its order through Theorem 14.121 We 
will show in the Appendix that .ft(A) is an indicator of the loss of accuracy in 
the computed solution of the linear system associated with FEM, GFEM, and 
SGFEM. ■ 

5 Applications: 

In this section we will present the SGFEM, when applied to three specific appli- 
cations. We will primarily address in detail the scaled condition number of the 
stiffness matrix of the method and show that the assumptions presented in the 
last section hold. The SGFEM, applied to each of these applications, will based 
on the uniform mesh {Tk}k^i\{o} with the set of vertices T, defined before. 
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5.1 Interface Problems 



Let a{x) in (|2.ip be a piecewise constant function and let / be smooth. We will 
consider two situations, namely, a{x) = ai{x) and a{x) = 02 (x), where 



ai{x) 



i, < a; < 5* 
1, b* <x<l 



and a2ix) = 



1, 0<x <bl 
i, hi <x <b^ 
1, 6^ < a: < 1 



We note that the solution u of (|2.ip does not belong to H^{^). 

We first consider a(a:) = ai(x). We consider the set 72 C T as before. There 
exists an m such that b* Gt,„_|_x= (xm, a;„j_(-i) and therefore, b* S H uj^n+i- 
For e T, we consider Vi = span{l, 1^9^*' = ^{\ / ai{t))dt\ . Clearly, for 
i 7^ m,m + 1, we have Vi = span{l,(a; — Xi-i)}. Therefore recalling that 
Vi — spanj^^j^"'}, where Ipi^ = (p'^i — S^t Vi = {0} for i 7^ m, m + 1. 

We set 72 = {xm, Xm+i} C T and from the definition of 52, we have 

^2 = ^ N,V, = 7V„F,„ + N„-,+iVra+l ■ 



We further note that (p ™ is linear on and therefore, y'l |t,„ ~ 0. Sim- 
ilarly, ^i""'^^' |t„,+2 ^ 0. Therefore t„j+i is the only enriched element, i.e., 
A^enr = {Tm+i}, and A22 = ^2™^"'^'. Also, wc Can easily show that ^i^'lr^+i = 
_[^m+i] 1^^^^^ ^ ^* _ x^+j3h with < /3 < 1. Then from a direct computation, 
we have 



A 



(m+l) 
22 



/i/32(l - /3)2(1 + 4/3)/6 - /3)(1 + 2/32)/3 



(5.1) 



Clearly, ^ 



(m+l) 



is associated with the vertices Xm, Xm+i- We choose the diag- 



onal matrix = diag{(5i'"'^'^ (5^""^'^}, where 



^(m+l) ^ ^-1/2^-1/2(1 - /3)-l, 4"+^' = ^ /3)-l/2. (5.2) 



Then 



_ £)(m+l)^(ni+l)£,(m+l) 



(|+/3-2/?2)/3 /3l/2(l„^)l/2(l + 4^)/6 



/3i/2(l-;3)i/2(i+4^)/6 



(l + 2/32)/3 



(5.3) 



The diagonal elements oi A'^^^^ arc 0(1) for all < /3 < 1. Also the eigenvalues 
of ^2?"''^'' are Ai = (2 — /3)/6 and A2 = (1 + /3)/2. Therefore, recalling Remark 
I4.14[ we have 

i||[^(m+l)]-l^||2 <^T^(™+1)^< ||[^(m+l)]-l^||2^ VxGR^, (5.4) 
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and hence, Assumption [3] is satisfied with L2 = | and U2 = I. 

— 1/2 

We set D2 = diag{o?i, ^2} with di ~ (A22)jj ■ Clearly, the diagonal ele- 
ments of A22 = D2A22D2 are equal to 1. Recall that 72 = {xrn.,Xm+i} and 
K-m = ^rn+i = {Tm+i}- Therefore, l{m,m + 1) = 1 and l{m + l,m + 1) = 2, 
where the index l{i,k) for Xi € T2 and S JCi was defined just before (|4.18[) . 
We also have 1C%^ = /C*j_|_i = {rm+i}- Therefore from (|4.18p . we have /S.„i = 
|-^(m+i)j_2 Am+i = [(Jj'"^^']"^. Also the vertices XmiXm+i G T2 arc as- 
sociated with the diagonal elements (A22)j„j„, (A22)j,„+ij„+i , respectively, of 
A22, where jm = 1, Jm+i =2. It is easy to check that 

1 < (A22)n A„, (A22)22'A„+i < 6 

and the Assumption [3] is satisfied with L3 = 1 and C/3 = 6. 

We have shown in Remark 14.141 that the Assumption [1] is always satisfied 
in 1-d; in this case Li = ^ and Ui = 2. Therefore, from Theorem 14.121 we 
have that A{A) = 0{h~^), and thus the GFEM with S = Si +S2 is indeed 
an SGFEM. We further note that Assumptions [1] [21 [3] are satisfied for any 
< /3 < 1, i.e., the constants Li,Ui, L2,U2, L3 and C/3 are independent of /3. 
Therefore M.{A) = 0{h~'^) even when /3 w or /3 « 1, i.e., when the point of 
discontinuity h* of ai(x) is close to the one of the vertices Xi (see also Remark 

EH). 

We next consider the (|2.ip with a{x) = 02(2;). We again choose Vi = 
span{l,(^]'^ = ^{\ I a2{t))dt\ . If the points of discontinuity h\, of a2{x) 

are separated, e.g., b\ £ti and 63 Gtj* with |Z — > 2, then we can again show 
that the GFEM with 5 = 5i + ^2 is an SGFEM, based on the arguments given 
above. 

Suppose there is an m such that 6| eT„j and 62 S'''m+i- Moreover, suppose 
6^ — Xrn-i + h/2 and b'2 = Xm + Ph with < /? < 1. Note that h\ is away from 
the vertices, whereas, could be close to either Xm (/? ~ 0) or Xm+i {P ~ 1). 
As before, let Vi = span{^!j*'}; clearly, Vi ~ {0} for i 7^ m — 1, m, m -I- 1. 
Therefore T2 = {x^-i, x^, Xm+i} and 

i^m— l,m,m+l 

We further note that = ^i"^^^!-^ = 0- Also it can be shown 

that '(^^l^ — and — "^i^^"^"*^^ |^ Therefore fCenr = 

} (i.e., T„i,Tm+i are the only enriched elements), and hence A22 is 
assembled from local stiffness matrices A^"^ and Al^^^^ . 
From direct computation, we get 

A A " 

16 32 

A A ' 

32 32 - 

and it is associated with the vertices Xm-\ and Xm- The matrix ^jT"*"^^ 
same as in (|5.ip and is associated with Xm and Xm+i- We choose Z?*^™) = 



A, 



(m) 
22 
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diagi6["'\ 4"^) with 



with 5["'^^\ Si""^"^^ as given in Then 



Si""^ = /1-1/2 and = diag(4™+'\ 



(m+2)N 



1 
32 



J_ 

f 

16 



Clearly the diagonal elements of ^2™'' are 0(1). The eigenvalues of ^2™^ are 
Ai = 1/32 and A2 = 3/32 and therefore (recall Remark l4.14p . 



32 



(m)l 



32 



(m)l 



(5.5) 



Next, the matrix Aj™"''^'' 



^(m+i)^(™+i)^(m+i) is same as the matrix 



(m+l) 
22 



given in (|5.3p . The diagonal elements of A 
are Ai = (2 - /3)/6 and A2 = (1 + /3)/2. Therefore 

1 



are 0(1) and its eigenvalues 



6 



(m+l)l-l. 



|2 <x^4!;^+^' 



X < ||[£)(™+l)]-lxj|2^ VxeR2. 



Thus the above inequality together with (j5.5p implies that Assumption [2] is 
satisfied with L2 ~ and U2 ~ for all < /3 < 1. 

The matrix A22 is assembled from the matrices Aj™^^'' and is given 

by 



k-22 



A 

16 

A 

16 





A 

32 



A 

32 



ft/3(l~/3)" 
3 



(§ + /?- 2/32) 
■(i+2/3) 



(1+2/3) 
M!^(i + 2/32) 



We choose D2 = diag((ii, (i2, ^3) with — (A22)i,j • Clearly the diagonal 
elements of A22 := D2A22D2 are equal to 1. Consider the vertex Xm & 72- 
Then /C,„ = {"Tm, Tm+i} and lim^m) = 2, l{m,m + 1) = 1. Also in this ease, 



JC* ~ JCi- Therefore, from (j4.18p . we have 
we can show that A„j_i = [(jj'"'']"^ and A 



m+l 



[6, 



(m+l)i 



(m+l)i_2 



Similarly, 



We also note that 



the vertices a;,„_i, x^, a;,„+i S T2 arc associated with the diagonal elements 



(A22) 



J, respectively, of A22, where jm-i 



(A22)j,„_ij,„_i, (A22)j„j„, V^22;j,„+ij„+ 

I7 Jm = 2, jjn+1 ~ 3. An easy calculation yields 

1 < (A22)n A,„_i, (A22)^2^A,„, (A22)^3^A,„+i < 16. 

Thus Assumption |3] is satisfied with — 1 and C/3 = 16 for all < /? < 1. We 
have shown before that Assumption [1] is always satisfied in 1-d. Therefore from 
Theorem l4.12l we infer that A{A.) = 0{h~'^); the result is true even when /3 « 
or ^ w 1. Thus the GFEM with 5 = 5i + ^2 is indeed an SGFEM. 

We remark that for a{x) = ai{x) or a{x) ~ 02 (x), we can show that there 
exists € Vi such that \\u — T^iU — ClUcc^i) = 0{h) for each Xi € 72- Thus 
using the standard interpolation estimates and using Theorem 14.71 we have 
\\u — Uh\\£{n) = 0(h), where Uh is the SGFEM solution. 
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Remark 5.1 Note that A'^'^^'^ and thus A22, A degenerate as /3 — >■ or /3 — >■ 1 . 
Let eo be small, say, eq = 10""'^^. We adjust the implementation when (3 < eo 
or 1 — ^ < eo by setting f3 ^ eq or I — f3 — eo, respectively. We emphasize that 
^(A) is bounded independently of f3. ■ 

5.2 Problems with singular solutions 

Let a{x) = 1 in ()2.ip and suppose f{x) be such that the solution u of p.ip - 
(|2.2p is of the form u = x" + uq, where ^<a<|,a^l, and uq is smooth 
with uo{0) 0. Clearly u ^ H'^{n). Let < £> < 1 and set (0,1?), 
fir ■■= {D,l). Then u € H^{nr) and |u|ff2(o^) < C[ |a;"|^/2(n^) + \uo\m{nJ- 
Clearly, |m|h2(q^) depends on D and is extremely large for Z? 0. 

We consider 7i C T as before. Let T2 ■— {xi e T : n il; 7^ 0}, where the 
patches uji have been defined before. Clearly, xo,xi G 72. Let A:* G / be the 
largest index such that g 72 for < i < fc* — 1. For € 72, let 

Vi = span{l, (ySi' = (x - Xj), 1^2' = 

and for Xi G T\T2, let Vi = span{l,(^^'' = (x — x^)}. Clearly, = {0} for 
Xi S T\72- For Xi e 72, we have 

= span{^^*' = crj 7^ 0, 

where (t,; := (x" — T^.x"')\i^r^ recall Xi^^x" is the piecewise linear polynomial 
that interpolates x" at the vertices {xi_i, Xi, x^+i} of for z 7^ 0, and IuiqX" 
interpolates x" at {xo,xi}. For an element C Zoi (with Xi € 72), we define 
a^^^ := (x" — /A;x")|rfc, where /fcx" G V^{Tk) interpolates x" at Xfe_i, Xfe. 
Clearly, X^.x" = /fex" on Tfc C aJ^. It is also clear that [a^'^^]' ^ on C uJi- 

We define ^2 = Z^iLo ^ -^'^jl^i and we consider the GFEM based S = Si + 82- 
We first address the convergence of the GFEM solution Uh- It is easy to show 
that for X, G 72, there exists ^, G Vi such that 

Also for Xi G T\T2, from standard interpolation result we have 
Therefore, from Theorem 14. 7[ there exists v Cz Si + S2 such that 

<C/i[K|1^2(0) + K|l,2(0„)]'^'. 

Thus we have — Uh\\£(n) < Ch, where Uh is the GFEM solution; note that C 
depends on |x"|/f2(-f2^-) and thus on D. 
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We note that f2; is independent of h. However, if _D = /i'*', 7 < 1 (i.e., 
= h'^), then one can show that ||u — Uh\\s{n) = 0{h^^^). Thus if we enrich 
only a fixed number of patches in the neighborhood of the singularity, we loose 
the optimal order of convergence. 

We now address the scaled condition number of the stiffness matrix A of 
the GFEM. We note that the matrix A22 is assembled from element stiffness 
matrices for the element r^, where Tk is enriched. We note that the set of 
enriched elements is given by ICenr '■= {tu € {ti\i£I\{o} '■ G 72}- We further 
note that if G /Cenr, then Tj € ICenr for ^ < j < k. Also from the definition of 
fc*, it is clear that t^-. S ICenr and Tj ^ ICenr for j > fc* + 1. Now, for G ICenr j 
k 7^ fc*, the matrices A*"^ are of the form A''^^ ~ {^hn}f m=i'^ entries b['2 
are as given by 



C= / {Nk-2+l<jy {Nk^2+,n<jy dx . 



Also, since ccfc. ^ T2, we have A^2 "* = [^11 "*] (an 1x1 matrix), where b[\ ' is 
given by the above expression. 

(k) 

Lemma 5.2 The entries of the matrix A22 are as follows: 



JTk JTk 

^^2^ =4? = / Nu-iNka'^dx. 



The proof is easy and we do not present it here. □ 

It is clear from above that for Tfc e ICenr and k ^ k* , the diagonal elements 
h'^ii ^ h^22 ^ ^ ^'^'^ therefore ^22' associated with x^-i, Xk- Also b[\^ > and 

(k*) 

thus A22 is associated with Xk*-i- A simple observation yields that the size 
of A22 is fc* X fc*. 

Let Tk G ICenr and set Xk-1/2 ■= {k — \)h\ Xk-1/2 is the mid-point of r^. We 
define 

Gk = \[xy{xk-i,2)\ = Ha - i)(fc - \r-^h^-\ 

Note that for 1 < j < fc* — 1, Tj+i G ICenr implies Tj € ICenr, and we have 



1<^ = ^ <3^-". (5.6) 

We now obtain a few results, which will be used to establish that the GFEM 
with 5 = 5i + 52 is an SGFEM. 

Lemma 5.3 For Xk G ICenr, there exist positive constants Cl,C2, independent 
of fc and h hut may depend on a, such that 

Gk ~ ^ 
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Proof: (a) First let 2 < fc < fc* and let g{x) = [cr('=)]'(a;) for x G Tfc. Then 
niax|5'(a;)| = max\[a'-''^]" {x)\ ^ \a{a ~ l)\x'^Zi 



= |a(a-l)|(fc-ir-2/."-2(^) 
= G.(^) <G.(-) :=M.. (5.7) 

Similarly, 

/k — -\ 2-ct /3\ 2-a 

mm\g'ix)\=Gk[^^) >Gfc(-) := m^. (5.8) 

We next note that g' does not change sign in Tk and thus g is monotonia in r^. 
Also, since /^^ gdx = a'^^^ \^ = 0, there exists a unique a;^ = Xk-\ + C,h ^ 
with < C < 1 such that g{x%) = and is characterized by 

g\dx^ I \g\dx. (5.9) 

We now obtain bounds on C, independent of fc. Since g^x^,) = 0, it is clear from 
the mean value theorem, (j5.7p . and (|5.8p that 

mk\x-xl\ < mm\g'\\x-xl\<\g{x)\ 

< ina:ii\g'\\x ~ xl\ < Mk\x - xl\, x e Tk . (5.10) 

Consequently, 



mfc^< / \g\dx<Mk^ (5.11) 
and 

Now from ((g?TT|) . and ([ET^ . we have 

mk ^ < / \g\dx= \g\ dx < Mk- ^' 

and thus 

VMfe _ (3/2)(2-")/2 



2 



Mfc + (3/2)(2-")/2 + (3/4)(2-a)/2 ' 

where we used the definition of and Mk given in (j5.8p and (j5.7p respectively. 

Using a similar argument we obtain a lower bound of C; we summarize the 
bounds of ( as 

< C < Cr, where 

(3/4)(2-")/2 



(3/2)(2-")/2 + (3/4)(2-a)/2 



(3/2)(2-")/2 

^"""^ " (3/2)(2-")/2 + (3/4)(2-")/2 ■ (^-^^^ 
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Finally, from (|5.10p and using the definition of (given in (IS-Zp ). we get 
/ \g\^dx = / \g\^dx+ / \g\^dx 

J Tk JXk-l Xl 

r2 



< / {xl - xY dx + Mi / {x - xlY dx 

Jxk-i J xl 

and similarly, we have 

\9\'dx > !!!^[c^ + [1 - cY] 

G| (3/4)2(2-") 

- 3 [Cf + (1-Crf]. 

Thus 



^.^3/2< HM^^rilL-(r.) ^^.^3/2^ for2<^<iV, 



where Ci* = (3/4)2- VICf + (1 ^ Cr-)3]/3 and Q ^ (3/2)2- V[C^ + (1 - CO^^IA 
(b) We now consider k = 1. We note that on ti = (0, /i), 

[cr(i)]'(a;) = ax°'-^ - h°'-\ 

By a direct computation, we get 



:= C* 







2a - 1 
Therefore, 

/J'|[(T(l)]f _ (Q-l)2fe2"-l _ 



Gl (2a- l)a2(a- l)2;i2a-424-2a (2a - 1) a2 2^-2" 



Hence we get the desired result with = min(Gi , G*) and Gj = max(G2 , G*). 
□ 

Lemma 5.4 Suppose G /Cenr I'^c^ l^t lk{x) he a linear function, defined on Tk, 
such that lk{xk-i) ~ yi and lk{xk) = Ui- Then there exists a positive constant 
G|, independent of k and h but may depend on a, such that 

\\[a^''^]'lk\\L^(r,)>ClGkh^/\yl+ylY'''- 
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Proof: (a) Let 2 < k < k* and define g{x) = [a^''^]'{x) for x S t^. On t^, we 
fiave seen in tlie proof of Lemma 15.31 tliat g{x\) = where x^, = x^-i + 3'iid 
0<Ci<C<Cr<l- We have also seen that 

mk\x-xl\<\g{x)\<Mk\x-xl\, V a; e , (5-14) 

where 

nik - Gfc(3/4)2-", Mk - Gfc(3/2)2-". 

Let 



1 - Cr 

Xfc = Xfc-i + Cr-h H — h. 



Then 



N/c |a;fe_i +Cr/i+ ^— ^^~2;fe_i - (ft-l > ^—^h. (5.15) 

Also from the definition of x^, it is clear that g(x) 7^ in {x^, x^) and thus from 
and dnHni), we have 

> mfc|xfc - x^l > — ^nikh. (5.16) 

Therefore, 

\g\'\hfdx> \gf\lkfdx>mlh' ^ / l/fepda;. (5.17) 

We make the change of variable y — ^-^jz^- Then 

where 

Hy) = ?fc(a;/c + — ^ — ) = 2/1 —^(1 - v) + y2 (^— + -^^y)- 

Thus F{yi,y2) is independent of fc. We next note that i^(?/i, 2/2) is a continuous 
function and F(^?/i,/3y2) = F{yi,y2). It is well known that the minimum of 
F(yi, 2/2) is attained on the compact set + yf = 1. Hence there is a constant 
Cmin, independent of k but may depend on Cr, such that 

il-Oh pMZPdx 
< Cl^J^^ < F{y,,y,) = . 

^ yi + i/2 

Thus from (|5.17p . we have 

" \9?\h?dx > Cl,^mlh'^l-^{yl+yl) 
= BfGlh^{yl+yl), 
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where Bf = (3/4)2(2-)c^„^(l - Cr?/8. 

(b) We now consider fc = 1. On ri = (0, h), we have g{x) = ax" ^ — h" ^. 
It is easy to see that g{x'^) ~ 0, where x\ — Qi with C = (^{o) ~ a^l^""-. 
Since (^{a) is increasing for ^ < a < | (with C redefined for a = 1), we have 
C < C ^ C(3/2) = (2/3)2. 

Set x\ — C*h + (1 — (^*)h/2. Since \g{x)\ is increasing in (xl,h), we have 
= |g(a;i)| on {xi,h). Therefore, 

'\g\''\h\^dx> r\g\'\li\'dx>gl,^ t\h?dx 



xi 



a2|a- l|2/l2("-2)/22(a-2) 

where 

o4-2qt:2 o4-2ar (1+C*) 

(j2 _ ^ gmm _ ^ I" 2 



a2(a- l)2/i2(a-l) a2(Q,_l)2 

As before, we can show that 



\h\''dx>C^J^^{yl + ylY/\ 



and therefore, 



r\9?\h?dx>BfG\h\yl+ylf/\ 
Jo 

where = C^Cmini^ — Q*)h/2. Finahy, defining C3 = min(_BJ', iJj) and 
recalhng that g = [c^*''']', we get the desired resuh. □ 

Now, for k ^ k* , consider the diagonal matrix D^'''' ~ dia.g{S[''\ 62''^ ) with 
(5f ^ = = Gj:i/i-3/2 and set i^2^ = £)('=)A^2^D('=). The diagonal elements of 
A22'' (see Lemma [5?^ are 

7(fe) _ /" ^2 ^'2 , _ 1 f A^2 /2 , 

Using Lemmas 15.41 and 15.31 it is clear that 

C;<P^^<-^I^ a'^dx<Cl 
where C2 , C3 are independent of k and h. Similarly, 



Cl<b^^<^j^ a''dx<C;. 

We let D^*^*) ~ [6['' ^] with S['^ ^ — G].}h^'^^^. Using similar arguments we 
show the CI < hfp < C2*, where i^f = D^''") Ai^P D^'''^ = [&<f Thus the 
diagonal elements of ^22"* are 0{1) for all € /Cenr- 

We next show that the element matrices A22'' satisfy the Assumption [21 
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Proposition 5.5 For Tk € fCenr, the matrices satisfies Assumption\^ 

Proof: Suppose k ^ k* and let x ~ {xi, X2)^ S M^. Then, using Lcnima l5.2[ we 
have 

x^A^2^x = b['\'xl + 2b[''^xiX2 + bi''^xl 

[xlN^_, + 2xiX2Nk^iNk + xlN^][a^''Y dx 

= 1 [.Ti7Vfc_i+X2iVfe]2[fTW]''rfa;<2(x?+x2) / 



and using Lemma 15.31 we have 



Next from Lemma 15.41 it is immediate that 

> C*Glh'{yf + yl)^C:\\[D^'^]-'^\f. 

Similar bounds fox ^x'^ for all a; G K also hold. Thus Assumption [2] is 
satisfied with L2 = Ci and C/2 = C2. □ 

Next, recalling that A22 is k* x k* , we choose the diagonal matrix D2 = 
diag(di , (^2 , • • • ,dk*) with dj ~ {A.22)Jj'^^ ■ Clearly, the diagonal elements of 
A22 — D2A22D2 are equal to 1. Note that (A22)jj, 1 < i < fc*, is associated 
with the vertex xj^i G 72- Also, (A22)ii = bfi and (A22)jj = ^22^^^ + for 
2<j< k*. 

Now, for Xi E T2 and i ^ 0, fc* — 1 (recall x.^ ^ T2 for k* < i < N), we 
have K-i = {Ti,Ti^i}, l{i,i) = 2, l{i,i + 1) = 1 and IC* = /Cj. Also K-q = {ti}, 
^(0,1) = 1, /eg = /Co and ICk'-i = {r^.}, l{k* - l,fc*) = 1, /C^._j = JCk'-i- 
Therefore, from (l4J8l) . we have = [4*^]"^ + for Xi G 7^ and 

i 7^ 1, fc* - 1; also Ao = [^^V^ and /^k'-i ^ [-5^ 'V^- 
We now show that Assumption [3] is satisfied. 

Proposition 5.6 Let {A.22) jj, 1 < i < fc*, be the diagonal elements 0/A22 and 
consider A^ for Xi G T2, defined above. Then Assumption\^is satisfied. 

Proof: Let Xi £ T2 and i 7^ 0,fc* — 1; is associated with (A22)jiji, where 
ji = i + 1. Therefore using the definition of (A22)i+i.i+i, d[^~^^\ and (52*\ we 
have 



(A22)- A. = (A22)- ,.,,A,: = + . (5.18) 
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Now using Lemmas 15.41 15. 3[ and (|5.6p , it is immediate that 



^ "22 ^ ^*2 



'-'a ^ "11 "-^i+i _ "11 



32(2-.) - G2^^;,3 Gf G\^,h^ 

and therefore, 



*2 

2 - "-'2 ! 



1 G2/l3 32(2"") 

< TTTT^TTTTTT < 7;::27^-;^7^ ■ (5-19) 



" 42^+61'^'' " Cf(l + 32(2-")) 
Similarly, we get 

1 Gl+ih^ 1 



Cf (1 + 32(2-)) - 6W+6(*^i) - 2C 



2 ■ 

2 



and combining (|5.18p . ()5.19p . we infer that there exist constants L3, t/3, such 
that 

U < (A22),^j\A, < C/3, 

where 



C/3 



26*2*^ C2*^(l + 32(2-")) ' 

1 32(2-a) 
2C3*2 ^ C3*2(l+32(2-")) ■ 



Thus Assumption [3] hold for Xi E T2, i Q,k* — 1. The proofs for Xi G T2, 
i = 0, fc* — 1 are simpler and we do not include them here, □ 

Based on Propositions 15.51 15.61 it is clear that Assumptions [2] and [3] are 
satisfied. Assumption [T] always hold in 1-d. Thus from Theorem 14.121 R{A) = 
0{h-^) and the GFEM with 5 = 5i + ^2 is an SFEM. 

5.3 Problems with discontinuous solutions 

We now address a problem, which is different from p.ip . Let = (0, 1) and set 
ri; = (0,c) and fir = (c, 1), where < c < 1 is fixed. Consider 

H{n) := |u e L2{n) : v{0) = v{l) ^ 0, j v'^dx < 00, 

and / v''^dx < 00 



Then {H{il), \\ ■ \\h) is a Hilbert space, where 
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We note that i?o(ri) C H{n) and functions in H{n) may have jump disconti- 
nuity at a; = c. 

For / G L2{^), wc consider the problem 

u e H{n), B{u, v) = F{v), W ve H{n), (5.20) 

where 

B{u,v) := / u'v' dx + / u'l/ dx a.nd F{v) :~ / fvdx. 
Jfii Ja,. Jn 

The biUnear form B{-,-) is coercive and bounded in H{il). Also F{-) is a 
bounded linear functional on H{n). Thus the problem (|5.20p has a unique 
solution. 

If / and the solution u e H{^1) of (|5.20p are smooth in 0/ and il,., then u is 
the solution of the boundary value problem 

— u" = / on ri(, — u" = / on 57^, 
m(0) = = and u'(c") = u'(c+) = 0. 

This problem mimics the problem with a crack in higher dimensions, where the 
solution is discontinuous along the crack line away from the crack-tip. 

We now give a characterization of the solution of (|5.20p . We will use the 
Heaviside function 

Lemma 5.7 Suppose u £ H{il) such that u'{c^) ~ u'(c+) = and 
{u"Ydx < oo, / {u"Ydx < oo. 



Then 

u{x) ~ s{x) + u{x), 
where s is a step function with discontinuity at x = c and u G H^{il) 

Proof: We first note that m(c~) and u(c"*") are well defined. We define 

u^u- ^ ^ ' [H,-l], (5.22) 

where Hc{x) is given in (|5.2ip . It is easy to check that 
"In, ="lo,' m(c-) =u(c-), 



Similarly, 



u\ 



m' „ , and v! {c ) — u' [c ). 



+ [u{c ) — w(c"'")], u{c'^) = u{c ), 



u'\^ =u'\^-^ , and u'{c )~u'{c ) 
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Note that u{c~) = u{c'^) = u{c~). We define u(c) = u{c~) so that u is contin- 
uous at x = c and thus is continuous on Q,. Also since w'(c^) = w'(c+) = 0, it 
is clear from above that u'{c~) = m'(c+) ~ 0. We define u'{c) = so that u' is 
continuous at x ~ c, and consequently u' is continuous in fl. Moreover, 



{u"fdx^ (u" fdx+ {u"fdx= {u" fdx+ {u"fdx <oo. 



n 

Thus u e H'^ifl) and considering s ^ _ i] ([s^J]), we get the 

desired result. □ 

Suppose c ^ T, i.e., c is not a vertex of the mesh. Therefore, there exists an 
m such that c S Tm+i and hence, c € Wm H Wm+i. We assume that m ^ ^, N; 
this is always achieved for h small. Since u{0) = u(l) = 0, we consider 71 = 
T\{xo,xn} (see Remark . 

For 1 < i < A^— 1, we consider Vi = span{l, (pi^ = (x—Xi), (^3' = Hc{x)} and 
we set Vq = span{(p^"' = (x — xq)} and Vn = span{i^^^' ~ {x — xn)}. Note that 
Vi G -ff(wi) for i e / (i.e., for x^ e T). Clearly, Vi = {0} for i e /\{m,m+ 1}. 
We set 72 = {a:^m, a;m+i} C T and define 

We consider the GFEM with S = Si+S2. 

Since (^2™' — He is constant in r™, we have ^2™' =0. Similarly, ^2"^^' lTm+2 
0. Therefore /Ccnr = {''m+i}- Moreover, the functions (^2™^ ^2™^^' ^'''^ discon- 
tinuous at a; = c, their values are zero at a; = Xm,Xm+i, and ^2™'!!- +1 = 

We assume that / is such that solution u £ H{n) of (|5.20p satisfies the 
assumptions of Lemma 15.71 and u = s + u, where s is a step-function with a 
discontinuity at a; = c and u g H^(fl). We now address the convergence of 
the GFEM solution. We first note that Theorem l4Jl hold for u € H{n) with 
£{Q),£{uji) replaced by H{D,), H{uJi) and with Vi e H{uji). Now for Xi € T\T2, 
we have u G H'^{uji) and from the standard interpolation result 



\u- 



-IloM\h(u,) = \\u-IujM\sioJi) <Ch\u\H2(u:,) = Ch\u\H2(cu,}- (5.23) 



For Xm G 72 , it is easy to show that there exists S Vm such that u — I^^ u — 
^ = u-l^^^uonujm- Therefore, ||it-Itj^,M-f™||i2(^^^) < C\ujm\ \\u-Ii^^^u- 
and from the standard interpolation result, we have 

\\u-I^^^u- ^"^Wniui^,) = \\u-^u,^u\\hHuj,„) < Ch\u\H^uj„,)- (5.24) 
Similarly, there exists ^,„+i G Vm+i such that 

\\u - Icj„, + iU - ^„i+l\\H{ur„,+i) < C!h\u\H2(^^^^^y 

Therefore combining (|5.23p . (|5.24p . (|5.3p and using the Theorem 14. 71 with mod- 
ifications as mentioned above, we infer that there exists v € S ~ Si + S2 such 
that 

I!" - v\\H{n) < Ch\u\H^n)- 
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Therefore, \\u — Uh\\H{n) = 0{h), where is the GFEM solution. 

We next address the scaled condition number of the stiffness matrix of the 
GFEM. Since T,n+i is the only element in /Cenr, we have A22 = Aj™^^'. 

Let c = Xm + I3h with < /3 < 1. A direct computation yields that 



A, 



(m+l) 
22 



(4-9;3 + 6/3^)/3 -i + 2;3-2^^ 
-i + 2/3 -2/32 (1 - 3^ + 6/32)/3 



Thus A^^~^^^ is associated with vertices Xm and Xm+i- We choose the diagonal 
matrix = dmg{6^^+'\ slf'^'^ } with 4"+'^ - 4"+'^ = /1I/2/2. Then 



_ £)(m+l)^(™+l)£)(m+l) _ 



The diagonal elements of A'^^^'^ arc 0(1) for any < /3 < 1. The eigenvalues 
of i^^"^^^ are Ai = | - 2/3 + 2/3^ - T and A2 = | - 2/3 + 2/3^ + T, where 
T = - Sip + 228/32 - 288/^3 + 144/34 (obtained from MAPLE). It can be 

shown that 

1 , 5 VU 1 , 5 

- < Ai < — , - < A2 < - + -^^ — . 

6- 6 6 2- 6 6 

Thus, as before, we have 



(4-9/3 + 6^2)/3 -1+2/3-2^2 
-i+ 2/3 -2/32 (1 - 3/3 + 6/32)/3 



1 ||[^(>n+l)]-l,||2 < x^4r'^X < + ^)||[i?(™+l)]-lx!|2, V X e M2, 

and the Assumption [2] is satisfied with L2 = ^ and [^2 = | + 

— 1/2 

We choose D2 = diag{di,d2} with di = (A22),j,; . Clearly, the diagonal 
elements of A22 = D2A22D2 are equal to 1. As in the first example (i.e., when 
a{x) = ai{x)) in Section \SA\ we have T2 = {xm,Xjn^i} and JCm = ^m+i = 
{Tm+i}- Therefore, l{m,m+l) = 1 and /(m+l, m + l) = 2. Also /C*„ = /C*j_|_i = 
{Tm+i} and therefore from (|4.18p . we have A™ = [(5i'"^^']^2 and Am+i = 
[^(m+l)j_2^ Also er2 are associated with (A22)y„j^„, , (A22)y„^.i2/„^.i , 

respectively, where ym = 1, and Hm+i =2. It is easy to check that 



3 24 

^ < (A22)ri'A™, (A22)22'A„,+1 < —. 

Thus Assumption [3] is satisfied with L3 = 3/4 and = 24/5. Therefore from 
Thcorcm l4.12[ we have .S(A) = 0{h~^), where A is the stiffness matrix, for all 
< /3 < 1. 



6 Conclusion 

The GFEM uses special enrichment functions, based on the available (or ex- 
tracted) information on the unknown solution of the underlying variational 
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problem. The use of special enrichment functions gives rise to the excellent 
convergence properties of the GFEM. In fact, for a given problem, it is possible 
to choose several classes of enrichment functions such that the GFEM, employ- 
ing each of these enrichment classes, will yield excellent convergence properties. 
However, GFEM employing some of these classes of enrichments could be ill- 
conditioned, i.e., there could be severe loss of accuracy in the computed solution 
of the linear system associated with the GFEM. The loss of accuracy could be 
much more than that experienced in a standard FEM. In this paper, we have pre- 
sented and analyzed a modification of the GFEM - the stable GFEM (SGFEM), 
which does not have the problem with severe loss of accuracy. SGFEM has all 
the advantages of the GFEM and is also very robust with respect to the parame- 
ters of the enrichments (e.g., the parameter /3 in Sections [5.11 and lO)) . The loss 
of accuracy is characterized by the scaled condition number and is expressed 
through Hypothesis H, which was validated based on various examples. 

The abstract framework developed in this paper has been applied to a one- 
dimensional problem for the clarity of exposition. This framework could also be 
applied to higher dimensional problems, which will be reported in a forthcoming 
paper. 
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7 Appendix 

Validation and Verification (V & V) is a fairly new field and is still in its devel- 
oping stage (pQ [71 [37]). Suppose a mathematical model of some "Reality" 
(e.g., a physical, chemical or biological system or process), formulated for a par- 
ticular goal or purpose, is given. The objective of V & V is to assess whether 
the predictions based on the computed solution of a mathematical model are 
reliable enough so that they could be the basis for certain decisions related to 
the goal. 

Validation is the process of building confidence on the mathematical model 
([Tl[321[37])- The process is of course is constrained by the cost, available time, 
and skills, as explicitly underlined in [35]. It is based on a set of properly se- 
lected problems and their mathematical models for which experimental data is 
available. These problems are called validation problems and they are chosen 
with varying level of complexity; more complex problems are closer to the "Re- 
ality" . Of course, obtaining the experimental data for the validation problems 
with increasing complexity is increasingly costly. The prediction based on the 
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computed solution of these problems is then compared with the experimental 
data. The assessment of the difference is based on a specified tolerance and a 
suitably selected metric (could be more than one) relative to the specific goal. 
If the measure of the difference is larger than the tolerance for any validation 
problem, the mathematical model is rejected. If none of the validation models 
are rejected, then one could have confidence that the mathematical model re- 
alistically describes the "Reality", with respect to the goal, beyond the scope 
of the chosen validation problems. The level of confidence will be based on the 
tolerance as well as the number and the selection of the validation problems. 
We mention that the set of the validation problems is finite, their selection has a 
large subjective component, and a philosophical question about the justification 
of the confidence in the mathematical model could certainly be raised (see [28]). 

Numerical algorithms and their properties obtained from the mathematical 
analysis are always based on various assumptions that are not satisfied when 
the algorithm is implemented on a computer. For example, infinite precision 
arithmetic is often assumed while describing a numerical algorithm or stating 
an inference about the algorithm obtained from the analysis. However, this as- 
sumption is always violated by a computer working with finite precision arith- 
metic. The output from the computer implementation of the algorithm may 
also depend, for example, on the package in which the algorithm have been im- 
plemented, the compiler, the processor, the computing platform with single or 
multiple processors, among other factors. Consequently, the output may vary 
even when the same outcome is predicted by the mathematical analysis for two 
different algorithms. For example, suppose the solution of the linear system 
Ax = 6 is sought using algorithms of the form P^APiZ = Pfb, x = PiZ, where 
Pi is a permutation matrix for i = 1,2. Both the algorithms should yield that 
solution X, however, the computed solutions could be different (see Problems 
la,b. below). Thus the implementation of a numerical algorithm in a computer 
is analogous to a "Reality" ; the goal is to obtain a particular quantity of inter- 
est for a particular purpose (related to a decision). The mathematical model 
of this "Reality" is the inference obtained from the mathematical analysis, or 
other statements based on the inference, about obtaining the quantity of inter- 
est from the algorithm. Therefore, the process of validation of the inference has 
to be performed to have confidence in the inference or a statement based on the 
inference. 

We have briefly formulated the following hypothesis in the Introduction. Let 
Ax = b, x,b G E" be a linear system, where the n x n matrix A belongs to a 
class of sparse matrices that include the stiffness matrices associated with FEM, 
GFEM, or SGFEM. Let x be the computed solution of the linear system, ob- 
tained from an elimination method, e.g., some variant of Gaussian elimination. 
Moreover, x is computed in finite precision arithmetic with machine precision e. 

— 1/2 

Let H = DAD where D is a diagonal matrix with Da = A^^ ' ; clearly. Ha = 1 . 
Recall the scaled condition number ^{A) of A is given by ^{A) := K2{H), where 
K2{H) = \\H\\2\\H~^\\2 is the condition number of H based on the || • ||2 vector 
norm. Also recall rj :— \\x — i;||2/||a;l|2- 
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Hypothesis H: For n, not small, 



T] w Cn^SiiA)e; /3 sa 0, 



(7.1) 



where x has been computed in an computing environment satisfying 
the IEEE standard for floating point arithmetic (with the guard 
digit), there is no overflow or underflow during the computation of 
X, and C, /3 do not depend on n as well as other factors mentioned 
before. 

The ~ in (|7.ip means there exist < Ci, C2 and < /? small, such that 
r] = Cnl^A{A)e with Ci < C < C2 and \/3\ < (3. Also this hypothesis addresses 
the range N for which not (almost) all digits of accuracy is lost (see Problem 
3a). 

Hypothesis H is based on certain mathematical inferences (results), which 
we will discuss later. The validation of (|7.ip with respect to the tolerance t = 
{'''i,T2} means that 6*2/6*1 < n and P < T2. Note that T2 is primary and should 
be small for confidence in (|7.ip . however ti could be allowed to be larger. The set 
of vahdation problems consists of stiffness matrices of FEM, GFEM, SGFEM, 
and other similar matrices, e.g., arising in finite difference method, applied to 
solve various linear elliptic variational problems of increasing complexity. For 
confidence in the Hypothesis, we require that (|7.ip is not rejected for any of the 
validation problems relative to the given tolerance t. We note that it is possible 
to select a tolerance such that the hypothesis is not rejected, however, the 
tolerance have to be admissible (e.g., reasonably small) for the decision making 
process. In our case, the decision will be whether to accept the SGFEM over the 
standard GFEM. We note that the class of matrices for which the hypothesis 
will be validated is not precisely defined, similar to a class of complex physical 
or engineering problem. 

Wc now give a theoretical rationale for (|7.1|) . There is a lot of literature 
available on the accuracy of the computed solutions of the linear system Ax = b. 
We particularly mention the classic [HI] and a modern book [55] with an excellent 
survey of the theoretical results in the area. Typically, the loss of accuracy in 
the numerical solution due to round-offs is analyzed by the backward error 
analysis. This analysis shows that the computed solution is the exact solution 
of a perturbed linear system, and it provides estimates of the perturbations in 
terms of the data of the linear system. A bound on the loss of accuracy in the 
computed solution, measured by rj defined before, is then obtained using the 
perturbation estimates. 

It is well known from a standard perturbation argument that for a full matrix 



where e is the machine precision and f{n) depends on the algorithm used to 
solve Ax = b (see e.g., [24l[27]). In Hypothesis H, we hypothesize that K2{A) is 
replaced by M.{A). We also hypothesize that Cin'^ < f{n) < €371^ and |/3| < ^, 
where C*i, C*2, and /3 are as defined before. It is important to note that in the 



A, 



V < f{n)K2{A)e, 



(7.2) 
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mathematical literature, only an upper bound of 77 is available; in contrast, the 
Hypothesis H addresses both the upper and lower bounds of rj. 

Consider the linear system DADz — Db, where D is a diagonal matrix with 
Da = 2^' in the rage of the floating point system. Clearly x — Dz. We now 
cite the following old result of F. L. Bauer ([5]): 

Theorem 7.1 Let x, z be the computed solutions of the linear systems Ax = b 
and DADz = Db, respectively, obtained from an elimination method with no 
pivoting. Furthermore, we assume that there is no overflow or underflow in the 
computation of x, z. Then all the digits of x andxjj are same, where xjj = Dz. 

We note that the result of the above Theorem is not true if the diagonal elements 
of I? are not binary. However in that situation, the quantities — — 
and |ji — are of the same order. 

We next note that it is possible to find a diagonal matrix such that H2{DAD) < 

K2(A). For example, for > 0, let A = ^ ^ . Then K2{A) = 1/^ 
is large for /i small. Let /.t = x'^ '^j where 1/2 < x < 2. Consider D = 



1 

2'*/2 



Then DAD = 



and 1 < K2{DAD) < 2 and conse- 



1 2'*/2 
X 

quently, K2{DAD) < K2{A) for /i small. Let D* be the diagonal matrix such 
that K2{D* AD*) ~ mini) K2(-DAZ?) (minimum over all diagonal matrices D 
with binary diagonal elements), then from the above theorem and (|7.2p . we 
have 

^ < f{n) K2{D*AD*)e < f{n) n2{A)e. 

Thus K2{D* AD*) provides more accurate information about 77 than K2{A). But 
in general, it is not easy to find either D* or H2{D*AD*). In Hypothesis H, 

— 1/2 

we used Si{A) ~ K2{DAD), where is a diagonal matrix with Da = A^^ 
and Da may not be binary. We note, however, that not using a binary only 
influences Ci, C2 by factors of 1/2 and 2 respectively. We also mention that 
in the literature ([Ml [2^), an upper bound of the form (|7.2p for rj is available 
with /(n) = Cn^ and K2{A) replaced by K{A) for symmetric positive definite 
linear systems solved by Cholesky decomposition. In Hypothesis H, we used 
f{n) = Cn^ , /3 « 0, based on our computational experience. 

We now consider a set of validation problems, whose exact solution (exper- 
imental data) is known. The solution to these problems will be computed on 
various computers using double precision, i.e., with 16 digits of accuracy. 

Problem 1: We consider approximating the solution u{x) ^ x oi the problem 
-u"{x) = 0, X e (0,l),it(0) = 0,u{l) = 1, by the FEM using piecewise 
linear finite elements. 

Problem la: We use the FE mesh vertices Xi = ih ioi i = 0,1, N 
and h ~ 1/N. The FE solution is same as the exact solution u 
of the problem. Let the associated linear systems be A'-'^'z'^^ = 
b^^\ The exact solution vector x*-^-* is known, namely, x^-^'' = ih, 



42 



z = 1, 2, • • • , iV. We will solve the linear system by the standard LU 
decomposition algorithm for sparse matrices without partial pivoting. 

Problem lb: We use the mesh vertices Xi = {N — i + l)h, i = 0, 1, • • • ,N . 
The FE solution is same as the exact solution u and let the associated 
linear system be A^'^'^x^'^^ = fe*-^^; it is known that x^P = {N — i + 
l)/i, i = 1, • • • ,iV. Note that the elements of x^^^ are the permuted 
elements of x'^^'> and thus ||a;(-'^^||2 ~ ||a;(^)||2. We will solve the linear 
system by the same algorithm as Problem la. 

The computations are performed on a Dell Latitude PC with INTEL 
C0RE(TM)2 CPU, 1.20GHZ. 

Problem 2: We approximate the solution u{x) = 1 of the problem —u"{x) = 
0, X G (0, 1), m(0) = 1, = 1, by the picccwise linear FEM based on the 
mesh vertices as in Problem la. Let the associated linear system be Ax = 
b. It is clear that the exact solution is given hy Xi = 1, i ~ 1,2, ■ ■ ■ ,N. 

Problem 2a: The linear system is solved by a sparse matrix direct solver 
superLU [29j on a single processor. 

Problem 2b: The linear system is solved by a sparse matrix direct solver 
MUMPS [2] on a single processor. 

Problem 2c: The linear system is solved by MUMPS, using parallel com- 
putation, on 128 processors. 

The computations were performed on the Lonestar system at Texas Ad- 
vanced Computing Center. Lonestar is a Linux based cluster comprised of 
1888 compute nodes connected via high speed quad-data rate infiniband, 
with each compute node containing two hex-core socket (INTEL Xeon 
5680 processors) for an aggregate system size of 22656 cores. Each core 
runs at a peak of 3.33GIIZ. 

Problem 3: We consider approximating the solution u{x) — x^ of the problem 
-v!'{x) ^ ~2, X € (0,1), u(0) = 0,u'(l) 2, by the GFEM based on 
S =^Si+S2 (see dS!])). We use Ui = 1 and ipf {x) = x^ , i ^ 0,1, ■ ■ ■ , N . 
We order the shape functions as NQLpf\ Ni, Niip^^^ N2,- ■ ■ , Njy, Nj^ip^^'^ 
and suppose the associated stiffness matrix is Ax ~ b, where A is of the 
order 2N + 1. The GFEM solution is same as the exact solution u. It is 
easy to see that X2i+i = 1, i = 0,2, ■, N and X2i = 0, i = 1, 2, • • • ,N. The 
linear system is solved by the same algorithm and on the same platform 
as in Problem la. 

Problem 4: We consider approximating the solution of the same problem in 
Problem 3 by the SGFEM based on S = Si+S2 (sec (|ITU1) ) with n, = 1, 
72 = T, and Tp^^^ = x^ —T^^.x^ . Wc order the shape functions as N^Tpf^ , Ni , 

Niip^i \ N2, - ■ ■ , Nn, A^Tv^i^' and suppose the associated stiffness matrix 
is Ax = b, where A is of the order 2N + 1. The GFEM solution is same 
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Figure 1: Log-log plots of 77''''' = ||a;''''' - :r'''' ||2/||a;'''' ||2 where i**^ is the computed 
solution of A'-'^'' x'-''^ = 6''"', k = 1,2, associated with FEM with vertices Xi = ih 
and Xi = {N — i)h, i = 0,1, ■ ■ ■ , N , respectively. 77'^^ rj^^^ have been computed and 
presented in (a) and (b), respectively, for A'' = 100, 200, ■ ■ ■ , 50000 



as the exact solution u and it is easy to see that X2i+i = 1, i = 0,2,-,N 
and X2i = {ihY, i = 1,2, ■ ■ ■ ,N. The linear system is solved by the same 
method and on the same platform as in Problem la. 

We will now validate Hypothesis H based on the validation problems de- 
scribed above. We will consider the tolerance r — {ti,T2), with ri = 400 and 
T2 = 0. 

Let and x'^^^ be the computed solutions of the linear systems A^^^x*^^^ = 
6'^-'^^ and A^^^x'^-' = b^^^ of Problem la and Problem lb, respectively. It can 
be shown that for large N, = ^{A'^^'>) « OAN'^. We have computed and 

presented the log-log plots of the relative errors rj^'^'^ = \\x^''^ — x*-'^-' ||2/||2;'''"'^ II2, 
k = l,2, with respect to = 100, 200, • • • , 50000 in FigurelH We have observed 
that ^ [0.4Ar2] < ^{k) < ^{k) jq 4^2] f^j. ^ ^ 1^2 with C2/C1 < 120 < n (note 
T2 =0). Thus we do not reject Hypothesis H. Note that we did not reject the 
hypothesis based only on the subset of meshes with the values of A^, mentioned 
above. Moreover, it is interesting to note that the plots of ry^^^ and ry'^^-' are quite 
different. Thus the computed solution is affected by changing the order of the 
FE mesh vertices, in spite of the fact that ||x*-^-'||2 = ||a;'^^^||2. 

In Problem 2, we solve the linear system Ax = b using two different soft- 
ware superLU and MUMPS; we also implement MUMPS on multiple proces- 
sors. Let be the computed solutions of Problems 2a, 2b, and 2c, 
respectively. These solutions were computed for 10 < < 10^, with 90 values 
of N in the range [10, 10^), with 400 values of in the range [10^, 10^), and 
360 values of TV in the range [10\ 10'+^), i = 3,4,5,6, and with A^ = 10^. We 
presented the log-log plots of rj^'^^ = \\x — i^*^^ ||2/||a;||2, k = a,b, for the values 
of N given before, in Figures [2^ and [2)d respectively. We observed that for 
A^ > 100, C2/C1 < 200 < Ti for both the problems. Thus we do not reject the 
Hypothesis H for Af > 100. Note that for N < 100, Figures EJi and b suggest 
that /3 « — 1 . It is also clear from Figure [^b that the implementation of the 
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io' 10= 10' 10* 10' 10' 10' 
(c) 



Figure 2: (a) Log-log plot of r;'"' = ||a:: — £ ||2/||a::||2 with respect to N, where 
a;'"' is the computed solution of Problem 2a (using superLU). (b) Log-log plot of 
j;'-''' = ||a:: — ||2/||2;||2 with respect to A'", where it") is the computed solution of 
Problem 2b (using MUMPS), (c) Semi-log plot of 100 * {rj^"'^ - r?^*"' )/??''^ with respect 
to A'", where 77''^' — \\x — 5:''^' ||2/||a;||2 and ft"' is the computed solution of Problem Ic 
(using MUMPS with 128 processors). Proportionally distributed 1931 values of N in 
the interval [10, 10^] are used in all the figures. 



algorithm in MUMPS changes drastically for N > 5 x 10'^; this is not the case 
with superL U, as seen in Figure [5^. Thus the computed solution depends on 
the software package, as mentioned before. For Problem 2c, we did not display 
the log-log plot of r;^'^) = ||a; — i('^'||2/||x||2 as it would be very similar to the plot 
of 77'-''^ in Figure [IJj. However, we computed Re = 100(77*^"^^ — j/''^)/?/'') - the 
"signed relative difference percent" — and presented the semi-log plot of Re in 
Figured}; for the same values of N, given before. For N < 5 x 10"^, we sec that 
i?e ~ and values of Re starts to oscillate for > 5 x 10'^. This indicates that 
the implementation in MUMPS changes drastically. Figure also suggests that 
rj^'^^ is larger than 77'''^ for most values on N, and 77'^^-' gets closer to t;^''^ as N 
increases. 

Let X be the solution of the linear system Ax = of Problem 3. We 
have M.{A) = 0{N^) (see Section 3.1). The log-log plot of 77 = - .t||2/||2:||2 
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(a) (b) 

Figure 3: Plots of ry = \\x — i||2/|l2;||2 where x is the computed solution of the linear 
system Ax — b of Problem 3, associated with the GFEM with vertices Xi = ih, 
i = 0, 1, ■ ■ • , TV. In (a), we used N = 50, 100, 150, • • • , 10000, and in (b), we used every 
value of A'' in the interval [100, 1000] to show the detail. 



with respect to = 50, 100, 150, • • • , 10000 have been presented in Figure [3^. 
In Figure ^jp, we show the details in the range 100 < N < 1000, where we 
have presented the log-log plot of 77 for every value of N in this range. Based 
on both these data (i.e., the values of 77 for every value of N in the range 
100 <_ iV < 1000 and for N = 1050, 1100, 1150, • • • , 10000),we have observed 
that CiTV* <r]< CaiV^ with C2/C1 < 340 < n (note T2 0). Thus we do not 
reject the Hypothesis H, again based on the subset of meshes with the values of 
N mentioned above. It is important to note that in Problem 3a, all the digits of 
accuracy were lost for N > 9000, and thus the Hypothesis H does not address 
the value of > 9000. We also computed 77 for every value of in the range 
9000 < N < 11000; 77 was of the order 1 and oscillated around 1. 

Let X be the computed solution of the linear system Ax = b of Problem 
4. We have shown in this paper that ^{A) = 0{h?). We have presented the 
log-log plot of 7; = \\x - x||2/||a:||2, with respect to A^ 50, 100, 150, • • • , 10000 
in Figured^, and for every value of A^ in the range 100 < N < 1000 in Figure 
[41d. Based on both these data (i.e., the values of 77 for every value of A^ in the 
r_ange 100 < A/_< 1000 andJV = 1050, 1100, 1150, • • • , 10000), we observed that 
CiN'^ < ri < C2N'^ with C2/C1 < 240 < n (note T2 = 0). Thus we do not 
reject the Hypothesis H (based on meshes with these values of A^). 

Thus we did not reject the Hypothesis H for any validation problems with 
respect to the tolerance ti = 400 and T2 = 0. But we would reject the Hypothesis 
H if we choose Ti = 300, since C2/C1 < 340 ^ ti in Problem 3. However, if 
the values of rj for every value of A^ in the range [100, 1000] were not available 
(see Figure [31d), then we will have C2/C1 < 250 < n, and we thus we would 
not reject Hypothesis H. Hence validation depends on the values of A^, i.e., on 
the number of validation problems considered, since each value of TV (in each 
of Problems 1, 2, 3, and 4) constitutes a separate validation problem. But as 
mentioned before, the choice of the tolerance depends on the type of decision 
related to the goal. For example in this paper, we have to decide whether to 
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, , 10^ 10' 10^ /|,\ 

(a) (b) 



Figure 4: Plots of ry = \\x — i||2/|l2;||2 where x is the computed solution of the linear 
system Ax = 6 of Problem 4, associated with the SGFEM with vertices Xi — ih, 
i = 0, 1, • ■ ■ ,N. In (a), we used iV = 50, 100, 150, • • • , 10000, and in (b), we used every 
value of A'^ in the interval [100, 1000] to show the detail. 



accept SGFEM over the standard GFEM. In this case, we may allow ri to be 
bigger; in fact if ti = 500, we still accept SGFEM over GFEM since the value 
of T] for GFEM wiU be much larger than the rj of SGFEM for large N. 
We summarize by stating that 

(a) we have confidence in Hypothesis H, based on the chosen val- 
idation problems (Problems 1-4). We underline that we have also 
considered other 2- and 3-dimensional validation problems for the 
Hypothesis H, which we do not present in this paper. We will present 
a more substantial validation of Hypothesis H in a future publica- 
tion. 

(b) Because of our confidence in Hypothesis H, we prefer the use of 
SGFEM over GFEM, since hnear system of SGFEM is less prone 
to the loss of accuracy than the linear system of the GFEM, when 
solved using an elimination method. 

Remark 7.2 As mentioned before, all the computations presented here were 
performed with 10~^^ accuracy. However, all the figures, presented above, indi- 
cate that that the apparent accuracy is about 10~^*. This is likely the effect of 
various cancelations. ■ 
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